重新审视Yliluoma的有序抖动算法

# 重新审视 Yliluoma 的有序抖动算法 来源：https://30fps.net/pages/revisiting-yliluoma-2/ **摘要：**本文讨论 Joel Yliluoma 在 2011 年提出的有序抖动算法（左图），比其他文章更详细地解释其内部原理，提出新的简化变体（中图），并将结果与最先进的算法（右图）进行比较。文章末尾附有源代码。 ## 1位简介 先打好基础。我假设你知道有序抖动是什么样子（如果不知道，请看上图）。在黑白图像中，编码很容易。首先，你需要从某个地方获取一个*阈值矩阵*。例如，一个 4x4 的拜耳矩阵，如下所示： \\begin\{bmatrix\} 0 & 8 & 2 & 10 \\\\ 12 & 4 & 14 & 6 \\\\ 3 & 11 & 1 & 9 \\\\ 15 & 7 & 13 & 5 \\end\{bmatrix\} 要对该矩阵应用于灰度图像，你需要遍历每个像素，找出它对应的矩阵元素（从概念上讲，矩阵平铺在图像上），并从矩阵中读取一个*阈值*。如果输入像素的灰度值高于阈值，则输出白色像素。代码如下： `` bayer_4x4 = np.array([ [ 0, 8, 2, 10], [12, 4, 14, 6], [ 3, 11, 1, 9], [15, 7, 13, 5]]) def dither_bayer4x4_1bit(img: np.ndarray): img_float = img / 255 # 在 [0,1] 范围内处理 H, W, _ = img.shape # 浮点输出图像 out = np.zeros((H, W, 1), dtype=float) for y in range(H): for x in range(W): color = img_float[y, x] threshold = (bayer_4x4[y%4,x%4] + 0.5) / 16.0 if color > threshold: out[y, x] = 1 else: out[y, x] = 0 return out `` 下面你可以看到它生成的结果。 对于 1 位输出，这种有序矩阵会产生相当粗糙的结果。如果让我选择，像 Atkinson 误差扩散抖动这样的方法在黑白图像中效果会更好。 但这只是热身练习。真正的问题是如何对彩色图像做同样的事情。为了创建类似这样的东西： 使用有序抖动的 16 色索引图像今天我们讨论的是如何精确地做到这一点。 ## 彩色有序抖动比你想象的既难又易 做好有序抖动出奇地棘手。但入门很容易。我们可以将阈值矩阵重新解释为结构噪声，将其平铺到整个图像上，并将其添加到原始颜色中，然后为每个像素找到最接近的颜色。1 位示例代码变成了这样： `` # 假设 "offset" 是两个调色板颜色之间的平均距离。 for y in range(H): for x in range(W): # 输入颜色 color = img_float[y, x] # (0, 1) 范围的阈值 threshold = (bayer_4x4[y%4,x%4] + 0.5) / 16.0 # 将阈值偏置到 (-0.5, 0.5) 范围，然后相加 moved = color + offset * (threshold - 0.5) # 找到最接近 'moved' 的调色板颜色 idx = find_index_of_closest(moved) inds[y, x] = idx `` 这种方法的效果实际上比预期要好，至少对于我的测试图像所使用的调色板来说是这样。唯一的问题是噪声强度，即上面提到的 `offset` 值。没有正确答案，但在我的实验中，这个神奇的公式效果不错： `` offset = dither_strength * 0.5 * median(pairwise_color_distances) `` 以下是这种“灰度偏移”方法的效果： 还不错！天空中的黑点有点粗糙，但可以通过降低抖动强度来清理。不幸的是，对于任意调色板，这种方法往往会产生去饱和的结果。但是，当调色板像这里一样是为图像设计的时候，效果还是可以的。 ### N 候选方法 解决*颜色选择问题*的一系列更高级的方法包括：为每个像素收集一组 *N* 个*候选*调色板颜色，为每个候选分配一个概率，然后通过随机方式或通过阈值矩阵选择一个。这些被称为 *N 候选*方法。我不打算在这里详细讨论，但会向你推荐一篇 2023 年的博客文章，标题是 *Ordered Dithering with Arbitrary or Irregular Colour Palettes* (https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/)，这篇文章解释了你需要知道的一切。请至少阅读其“概率矩阵”(https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/#n-candidate-methods)部分。 这些算法试图满足的一个关键属性是*局部均值再现*：当从远处看（或模糊）抖动结果时，它应该看起来像原始图像。这转化为最小化输入像素颜色 \\mathbf\{p\} 与所有选中候选的加权和 w\_1 \\mathbf\{r\}\_1 + \.\.\. + w\_N \\mathbf\{r\}\_N 之间的距离。更多细节，我建议阅读上面链接的博客文章。 ### Knoll 算法 16 色图像使用 Knoll 算法抖动。Thomas Knoll 的抖动算法非常好地最小化了上述距离，该算法在 Adobe Photoshop 中广为人知（毕竟他是其发明者）。Knoll 算法首先设置一个*目标颜色* \\mathbf\{x\}\_1 为最接近输入颜色 \\mathbf\{p\} 的调色板颜色。这个调色板颜色 \\mathbf\{r\}\_1 是第一个候选颜色。然后，算法测量 \\mathbf\{r\}\_1 和 \\mathbf\{p\} 之间的误差，并将目标向相反方向移动： \\mathbf\{x\}\_\{2\} = \\mathbf\{x\}\_1 + (\\mathbf\{p\} - \\mathbf\{r\}\_1)。 现在重复该过程，找到最接近目标 \\mathbf\{x\}\_2 的新颜色。这将是第二个候选颜色。它补偿了前一个候选的误差；如果找到的第一个调色板颜色太蓝，那么这一个将是偏黄的。经过 N 轮这样的迭代，算法访问了 \\mathbf\{p\} 周围的不同颜色。每个颜色可能被多次访问，最终候选概率与访问频率相关。该过程可以总结如下： Knoll 算法的误差补偿循环1. 找到最接近目标点 \\mathbf\{x\}\_i 的调色板颜色 \\mathbf\{r\}\_i。 2. 增加 \\mathbf\{r\}\_i 的权重，加 1。 3. 移动目标 \\mathbf\{x\}\_\{i\+1\} = \\mathbf\{x\}\_i + (\\mathbf\{p\} - \\mathbf\{r\}\_i)。 4. 重复 N 次。 5. 用 N 归一化权重。 当过程重复到 \\mathbf\{x\}\_\{i\+1\} 时，最近点查询将在 \\mathbf\{p\} 的另一侧进行。这样选择连续的点近似于一个凸包。我知道这个过程可能仍然显得抽象，但重点是算法既简单又有效。有一些微妙之处，比如按亮度排序 (https://matejlou.blog/2023/12/06/ordered-dithering-for-arbitrary-or-irregular-palettes/#appendix-i-candidate-sorting) 以保持像素间颜色选择的一致性，但在 Python 中归结为： `` # 在主循环外分配工作数组 weights = np.zeros(K, dtype=float) error = np.zeros(3, dtype=float) for y in range(H): for x in range(W): color = img_float[y, x] # 清空工作数组 weights[:] = 0.0 # 候选权重 weight_sum = 0.0 error[:] = 0.0 # 累积误差补偿 # 寻找候选，可能少于 N 个唯一颜色 for _ in range(N): idx = find_closest(color + error * dither_strength) weights[idx] += 1 weight_sum += 1 error += color - palette_float[idx] weights /= weight_sum threshold = (bayer_4x4[y%4,x%4] + 0.5) / 16.0 # 输出累积权重和首次超过从 4x4 矩阵读取的阈值的候选索引 # 如果调色板索引未被使用，其权重将为零 # 并且不会被此循环选中。 cumulative_sum = 0.0 for idx in luma_order: w_ko = weights[idx] cumulative_sum += w_ko if cumulative_sum > threshold: break inds[y, x] = idx `` 注意，抖动强度可以通过调制误差补偿的尺度来调整。你也可以研究 matejlou 的着色器版本 (https://www.shadertoy.com/view/dlcGzN)。这些背景知识应该足以理解接下来要介绍的替代方案。 ## Yliluoma 算法 [](https://30fps.net/pages/revisiting-yliluoma-2/yliluoma2-full.jpg)2011 年，Joel Yliluoma 提出了一系列抖动算法作为 Knoll 算法的替代方案。它们在他的网站上的一篇文章中进行了描述，标题为 *Arbitrary-palette positional dithering algorithm* (https://bisqwit.iki.fi/story/howto/dither/jy/)。据我所知，它们并没有获得太多关注。我仔细研究了它们，发现其中一个特别有趣。 ## Yliluoma-2 简化版 之前使用 Yliluoma 算法 #2 抖动的 16 色图像。我将专注于文章中介绍的*第二个算法* (https://bisqwit.iki.fi/story/howto/dither/jy/#YliluomaSOrderedDitheringAlgorithm%202)。确切地说，是其 C++ 实现中的变体。我们称它为 *Yliluoma-2*。 从高层来看，Yliluoma-2 与其他 N 候选算法（如 Knoll 算法）类似：在每一个像素上，从调色板中选择加权候选颜色，然后根据阈值矩阵（或噪声）输出一个。区别在于候选的选择方式，在这种情况下，是通过使用一个独特的误差公式测试每一个调色板颜色来完成的。 ### 不需要复杂的色差公式 有一种普遍的看法，认为高质量抖动需要复杂的色差计算。我认为这是错误的，你真正需要的是 (a) 对绿色通道赋予更多权重，以及 (b) 强调亮色之间的差异。 文章提出了一种亮度加权的色差公式（本网站以前 (https://30fps.net/pages/better-srgb-to-greyscale/) 讨论过），我实现了它。你可以在下面第一张图中看到它的结果。实现同样效果的一个更简单的方法是在抖动之前对图像*和*调色板进行去饱和处理。最终的索引图像仍然使用原始、多彩的调色板。例如，libimagequant (https://pngquant.org/lib/) 将 RGB 颜色乘以 (0.5, 1, 0.45) 并将每个通道提升到 0.8 次方。我在下面的第二张图中这样做了，效果同样好。 当然，如果你想要高质量的还原，比如“眯起眼睛看起来像原图”，那么应该在线性空间中进行抖动。我还想说，如果你的输出分辨率低，像素大，你可以做任何你想做的事。你正在构建一个看起来非常不同的东西，并希望观看者以与原始图像相同的方式解读它。 ### 通过点到线段最近点测试选择候选人 回到主要话题。我花了一个愉快的下午研究代码，并得出了它的几何解释。回想一下 N 候选方法是如何既找到候选调色板颜色又为它们分配权重（概率）的。Yliluoma-2 中的颜色选择循环会跟踪迄今为止选中的候选颜色的*指数移动平均 (EMA)* (https://en.wikipedia.org/wiki/Exponential_smoothing)。换句话说，在每次搜索迭代中，该例程通过 \\mathbf\{x\}\_\{i\+1\} = (1 - t)\\mathbf\{x\}\_i + t\\mathbf\{r\}\_i 用一个新的候选 \\mathbf\{r\}\_i 更新移动平均值 \\mathbf\{x\}\_i，其中 t \\in [0,1] 是一个每次迭代变化的混合权重。（也许 t 应该被称为 t\_i？） 好的，所以候选通过 **lerp()** 混合到运行平均值中。但是候选是如何选择的，混合因子 t 又是从哪里来的呢？候选点是通过测试每个调色板颜色 \\mathbf\{c\}\_k 与当前均值 \\mathbf\{x\}\_i 之间的线段来选择的。最接近输入颜色 \\mathbf\{p\} 的线段告诉我们选择哪个调色板颜色作为候选 \\mathbf\{r\}\_i。 t 的值是产生最接近输入的线段上的最近点的混合因子。这很难解释，但从视觉上看应该更容易理解： 从 \\mathbf{x}_i 到 \\mathbf{r}_i 的线段上最接近输入 \\mathbf{p} 的点由 t 标识。在这个例子中，t 大约为 0.6。总结起来：在所有调色板颜色中，选中的候选在与当前均值线性组合时最能逼近输入。 ### 指数移动平均循环的详细说明 下面是更详细的算法。理解它对于理解下面的内容并非必要，但我认为将其记录下来是明智的。 在第一次迭代中，候选被选为最接近输入颜色的调色板颜色：\\mathbf\{x\}\_1 = \\text\{closest\}\(\\mathbf\{p\}\)。其权重为 1。其余迭代过程如下： Yliluoma-2 的指数移动平均循环（简化版）1. 将每个调色板颜色 \\mathbf\{c\}\_k 都视为可能的候选。 - 找到最小化平方距离 ||((1-t)\\mathbf\{x\}\_i + t \\mathbf\{c\}\_k) - \\mathbf\{p\}||^2 的混合因子 t。 - 如果距离是目前找到的最短的，则： - 更新 t\_\\text\{best\} = t - 更新 \\mathbf\{c\}\_\\text\{best\} = \\mathbf\{c\}\_k。 - 更新 k\_\\text\{best\} = k。 2. 选择候选 \\mathbf\{r\}\_i = \\mathbf\{c\}\_\\text\{best\}。 3. 通过最佳混合因子将均值点移向候选：\\mathbf\{x\}\_\{i\+1\} = (1 - t\_\\text\{best\}) \\mathbf\{x\}\_i + t\_\\text\{best\} \\mathbf\{r\}\_\{i\} 4. 将 t\_\\text\{best\} 添加到候选 \\mathbf\{r\}\_i 的权重中。 5. 重复 N 次。 通过 `weights /= sum(weights)` 归一化权重。 与 Knoll 算法一样，Yliluoma-2 也最终近似于一个凸包。在每次迭代中，目标 \\mathbf\{x\}\_\{i\+1\} 被移动到之前找到的线段上的最近点。由于这个点位于旧目标 \\mathbf\{x\}\_i 和输入颜色 \\mathbf\{p\} 之间，它必须在 \\mathbf\{p\} 的对侧（笼统地说）。这就是目标（即运行平均值）围绕 \\mathbf\{p\} 移动的原因。 *在上面的解释中，我试图通过选择带有显式平均步骤的公式来使过程易于理解。原始 C++ 代码使用运行平均值，可能是为了速度。* ## 如何找到 *t* 上面，我没有说明如何精确找到线段上的最近点。最简单的方法是进行*参数扫描*：测试不同的 t \\in [0,1] 值，并保留最接近输入点的那个。下面是一个示意图。我们试图在线段 AB 上找到最近于点 C 的点 D。下面的红点代表我们正在测试的各个 t 值： 我们想要线段 AB 上最接近 C 的点。不幸的是，参数扫描不太可能命中精确的最优 *t* 值。它也需要时间，但如果你有一个复杂的色差公式，它可能是唯一可行的方法。但还记得我们是如何确定亮度加权可以通过预处理完成的吗？这就是为什么我们在搜索 t 时可以忽略感知因素，直接解决几何问题。 ### 解析解 *t* 假设 A \\neq B，我们得到 t\_\\text\{closest\} = \\frac{(C - A) \\cdot (B - A)}{||(B - A)||^2}, 其中 “\\cdot” 表示点积。由于我们想要 A 和 B 之间的点，我们必须将范围限制在 [0,1]： t = \\max(0, \\min(1, t\_\\text\{closest\})) 并且像之前一样，D = (1-t)A + tB。 代码实现一致`` c_k = palette_float[k] # B delta = c_k - mean # B - A t = 0.0 delta2 = delta.dot(delta) # ||B - A||^2 # 避免除零 if delta2 >= 1e-9: t = (p - mean).dot(delta) / delta2 # 将选择限制在给定范围内 t = max(min_t, min(max_t, t)) mix = (1-t) * mean + t * c_k # D diff = p - mix cand_dist = diff.dot(diff) # d^2 `` ### *t* 的限制 由于精确解可以位于任何位置……