image-rs 中 fast_blur 速度提升 5 倍

# image-rs 中的 fast_blur 提速 5 倍 – Arthur Pastel 来源：https://apas.tel/blog/optimizing-image-rs-blur 几周前，我为了找乐子，开始寻找可以优化的东西。Rust 的 `image` (https://crates.io/crates/image) crate 引起了我的注意：我之前用过它，它是下载量最高的 crate 之一，而且图像处理本身计算量就很大。我很快发现了一个叫做 `fast_blur` 的方法。光看这名字，就值得尝试优化一下。结果是：**在 `u8` 像素的图像上最多快 5.9 倍**。 fast_blur σ=50 ×5.1 更快（越低越好） x86_64 上的墙钟时间测量 在深入优化细节之前，我们先退一步，理解模糊是如何工作的，不同算法之间的权衡，以及 `fast_blur` 在模糊质量和性能谱系中的位置。 ## 模糊性能与质量 (https://apas.tel/blog/optimizing-image-rs-blur#blur-performance-and-quality) ### 高斯模糊 (https://apas.tel/blog/optimizing-image-rs-blur#gaussian-blur) 可能是最常见的模糊，也是人们提到“模糊”时通常想到的效果。它将每个像素替换为其邻域像素的加权平均，其中 σ 控制影响范围的大小。实际上，由于像素是离散的，这转化为一个 k×k 的卷积核，其中核大小 k 与 σ 直接相关。 原始图像 原始图像 高斯模糊，σ=15 sigma=15 的高斯模糊 由于卷积运算，朴素二维高斯模糊的每个像素复杂度为 O(k²)。但高斯模糊可以分解为两个一维 pass，从而将复杂度降至每个像素 O(k)。这正是 image-rs 的 `blur` (https://docs.rs/image/latest/image/enum.DynamicImage.html#method.blur) 实现的方式，但如果能牺牲一些质量，我们也可以用更快的速度实现模糊效果。 ### 方框模糊 (https://apas.tel/blog/optimizing-image-rs-blur#box-blur) 这种方法简单得多：每个像素被替换为其固定大小窗口内邻域像素的平均值。本质上，是一个具有均匀权重的卷积核。这确实能模糊图像，但内核缺少平滑性，产生的模糊效果比高斯模糊更生硬、更不自然： 方框模糊，σ=15 sigma=15 的方框模糊 高斯模糊，σ=15 sigma=15 的高斯模糊 与高斯模糊一样，方框模糊也可以分解为两个一维 pass。此外，由于所有权重都等于 1/k²（对于 k×k 的方框模糊），我们可以将除法提取出来：先累加所有值，最后统一除一次。更进一步，每个 pass 都可以使用滑动窗口，而不是从头计算平均值。最终，每个像素在每个 pass（垂直和水平）中大约只需 2 次加法，以及 1 次除法来平均求和，使得无论核大小如何，每个像素的复杂度为 O(1)。 因此，我们每个像素的复杂度从 O(k) 降到了 O(1)——快得多。代价是：结果比高斯模糊更不平滑、更不自然。 ### 快速近似高斯滤波（即 `fast_blur`） (https://apas.tel/blog/optimizing-image-rs-blur#fast-almost-gaussian-filtering-aka-fast_blur) 这是一种两全其美的方法，前提是不严格要求精确准确。其思想¹ (https://apas.tel/blog/optimizing-image-rs-blur#user-content-fn-1) 是：使用三个精心选择核大小的方框模糊依次应用，结果能紧密近似真正的高斯模糊。每增加一个方框模糊 pass，内核都会进一步向高斯平滑： 好消息是，实际上我们只需使用方框模糊，就能得到外观好得多的模糊效果。 方框模糊，σ=15 sigma=15 的方框模糊 快速模糊，σ=15 sigma=15 的快速模糊 高斯模糊，σ=15 sigma=15 的高斯模糊 由于这仅使用方框模糊，现在我们有了一个外观不错的模糊，每个像素的代价为 O(1)，与 sigma 无关，这正是 image-rs 中 `fast_blur` (https://docs.rs/image/latest/image/enum.DynamicImage.html#method.fast_blur) 的实现。 ## 热路径 (https://apas.tel/blog/optimizing-image-rs-blur#the-hot-path) 现在我们对算法、复杂度和质量有了更多了解，可以开始看代码并尝试让它更快。第一步是进行性能分析，看看时间花在哪里。我使用了我们在 CodSpeed (https://codspeed.io/) 构建的性能分析工具，以及 MCP 集成和智能体来快速识别热点。 `fast_blur` 在 `u8` 图像上的火焰图主要由三个开销主导： - **`roundf`：27%** 总时间 - **`to_f32`：22%** 总时间 - **`min`/`max`：20%** 总时间 热路径将每个 `u8` 子像素转换为 `f32`，用浮点算术累加，然后通过 `rounding_saturating_mul` 转换回来，后者内部调用了 C 库的 `roundf`。对于一张 1920×1080 的 RGBA 图像，有超过 800 万个像素，每个像素在 6 次方框模糊 pass（3 次水平 + 3 次垂直）中多次经历这种浮点数往返。 ## 将浮点运算替换为整数运算 (https://apas.tel/blog/optimizing-image-rs-blur#dropping-floats-for-integers) `u8` 像素的值范围是 0-255，一个方框模糊核最多累加 ~宽度（或高度）个像素。我们在 u32 中能存储的最大累加器值为 `u32::MAX / 255`，大约是 1680 万，相当于两整张 4K 图像。仅仅是模糊核就足够了。实际上，对于任何我们可能用到的模糊半径，这都有足够多的余量。所有累加都可以用整数算术完成，从而消除浮点转换、`roundf` 调用以及由 `rounding_saturating_mul` 引起的 `min`/`max`（后者用于将结果限制在 u8 范围内）。 去掉所有样板代码后，内循环从这样： ``` let mut sum: f32 = 0.0; for i in kernel_range { sum += src[i].to_f32().unwrap(); } *dst = rounding_saturating_mul(sum, 1.0 / kernel_size as f32); ``` 变成这样： ``` let mut sum: u32 = 0; for i in kernel_range { sum += src[i] as u32; } *dst = ((sum + kernel_size / 2) / kernel_size) as u8; ``` 这产生了完全相同的模糊输出，但无需经过昂贵的浮点运算。需要注意的是，除法之前添加了舍入偏置（`kernel_size / 2`）。在整数除法前加上这个值，能给出与正浮点数上 `roundf` 相同的舍入行为，但每个像素只增加一次加法。 ### 设计 trait (https://apas.tel/blog/optimizing-image-rs-blur#designing-the-trait) 该技术适用于 `u8` 像素，但 `image` 支持许多其他原始类型：`u16`、`f32`，甚至 `f64`。有些用户模糊 HDR 浮点图像，有些则模糊 8 位缩略图。强制对 `f32` 像素使用整数累加器毫无意义，反之，保持所有路径使用 `f32` 会阻碍我们刚刚为 `u8` 设计的优化。因此，优化需要与现有的浮点路径共存，并且不能在每个模糊 pass 中散布 `if pixel_type == u8` 分支，因为那样要么失去意义，要么使代码难以维护。 由于两种情况下操作的形式相同（从零开始，累加采样值，除以并存储），只是类型和算术不同。经过与维护者的几次迭代，设计最终落定为一个 `BlurAccumulator` trait：以源像素类型为参数，实现目标为累加器类型。去掉错误处理和边界情况方法后： ``` pub(crate) trait BlurAccumulator: Copy { type Weight: Copy; const ZERO: Self; fn from_primitive(value: T) -> Self; fn create_weight(kernel_size: usize) -> Self::Weight; fn to_store(self, weight: Self::Weight) -> T; } ``` 有了这个 trait，内循环对像素类型及其累加器都是泛型的： ``` let weight = Acc::create_weight(kernel_size); let mut sum = Acc::ZERO; for i in kernel_range { sum += Acc::from_primitive(src[i]); } *dst = sum.to_store(weight); ``` 通过将 trait 放在累加器上（而不是源类型），`u32` 知道如何接收 `u8` 采样并对它们求平均，而像素类型只是其中一个输入。两个实现覆盖了所有情况。`u32` 用于 `u8` 快速路径： ``` impl BlurAccumulator for u32 { type Weight = u32; // kernel size const ZERO: u32 = 0; fn from_primitive(v: u8) -> u32 { v as u32 } fn create_weight(ks: usize) -> u32 { ks as u32 } fn to_store(self, ks: u32) -> u8 { ((self + ks / 2) / ks) as u8 } } ``` `f32` 用于其他所有类型，保留原始浮点路径不变： ``` impl BlurAccumulator for f32 { type Weight = f32; const ZERO: f32 = 0.0; fn from_primitive(v: T) -> f32 { v.to_f32().unwrap() } fn create_weight(ks: usize) -> f32 { 1.0 / ks as f32 } fn to_store(self, w: f32) -> T { rounding_saturating_mul(self, w) } } ``` 源原始类型上的一个小型伴随 trait 在类型层面选择合适的累加器： ``` pub trait WithBlurAcc: Sized { type BlurAcc: BlurAccumulator; } impl WithBlurAcc for u8 { type BlurAcc = u32; } impl WithBlurAcc for u16 { type BlurAcc = f32; } // ... f32, f64, i8, i16, ... 都映射到 f32 ``` 模糊 pass 对两种类型都是泛型的，正确的实现在编译时完全确定。编译器为 `u8` 路径选择 `u32`，为其他所有类型选择 `f32`，并进行相应内联。实际上，由于当核大小超过 1600 万像素时 `u32` 累加器可能溢出，当半径极大时仍有一次运行时检查以回退到 `f32` 累加。但正常使用场景完全在 `u32` 限制之内，因此这不是问题。 ### 影响 (https://apas.tel/blog/optimizing-image-rs-blur#impact) 仅凭去掉浮点运算，所有 `fast_blur` 基准测试的墙钟时间就获得了 **~×1.83** 的加速，且测试套件中没有质量回归： fast_blur σ=50 ×1.8 更快（越低越好） 之前：f32 累加器（主分支）。之后：u32 累加器加整数除法。 ## 避免除法 (https://apas.tel/blog/optimizing-image-rs-blur#avoiding-divisions) 整数累加器优化替代了浮点运算，但留下了 `(sum + half_k) / kernel_size_u32`。这消除了 `roundf` 和 `to_f32`，但 `div`，即除法指令，仍然是一个显著的热点。即使我们讨论的只是一条汇编指令，整数除法也是现代 CPU 上最昂贵的算术运算之一。一条 `div` 指令会阻塞流水线 20–30 个周期² (https://apas.tel/blog/optimizing-image-rs-blur#user-content-fn-2)，并且与大多数算术指令不同，它不能被流水线化，意味着 CPU 会停顿直到完成。对于一张 1920×1080 的 RGBA 图像，经过 3 次方框模糊 pass，有超过 5000 万次除法伴随着停滞的流水线。但是，在单次模糊 pass 中，`kernel_size` 对所有像素都是常量，因此必然可能预先计算一部分工作。 ### 技术原理 (https://apas.tel/blog/optimizing-image-rs-blur#the-technique) 在实数算术中，同时乘除 2³² 不会改变结果： n/d = (n·2³² / d) / 2³² 只有当 2³² / d 恰好是整数时才有用，但这种情况很少。然而，Granlund & Montgomery (1994)³ (https://apas.tel/blog/optimizing-image-rs-blur#user-content-fn-3) 表明，使用倒数的上取整而不是精确倒数，得到的近似值保证向下取整到正确的整数商： ⌊n/d⌋ = ⌊ n·⌈2³²/d⌉ / 2³² ⌋ 只要 (d-1)·n < 2³²，这个结果就是精确的（在我们的使用中几乎总是成立）。由于我们总是知道 d（核大小）在编译时已知，并且 n 是有上限的（由最大可累加值决定），我们可以预先计算一次倒数，然后在循环中用它替换除法。这样做还有一个额外的好处：倒数乘法可以被流水线化得很好，而除法不能。 因此，我们用一条 `mul` 和一个移位替换一条 `div`，这是最简单的算术汇编指令之一，并且流水线化效果极佳。简化后，原来的代码： ``` *dst = ((sum + kernel_size / 2) / kernel_size) as u8; ``` 变成： ``` *dst = (((sum + kernel_size / 2) as u64 * reciprocal) >> 32) as u8; ``` 在每次模糊 pass 之前定义一次倒数： ``` let reciprocal = (u32::MAX / kernel_size) + 1; ``` 这将昂贵的 `div` 替换为一个 64 位乘法和一个位移，在现代 CPU 上大约只需要 3–4 个周期。 ### 实现 (https://apas.tel/blog/optimizing-image-rs-blur#implementation) 为了在 trait 内部清晰地实现这一点，`u32` 累加器的 `Weight` 同时携带倒数和舍入偏置： ``` #[derive(Clone, Copy)] struct U8Weight { reciprocal: u32, // ceil(2^32 / kernel_size) rounding_bias: u32, // kernel_size / 2 } ``` 并且 `u32` 上的 `BlurAccumulator` 实现变为： ``` impl BlurAccumulator for u32 { type Weight = U8Weight; // ... fn create_weight(kernel_size: usize) -> U8Weight { let ks = kernel_size as u32; U8Weight { // ceil(2^32 / ks) = floor((2^32 - 1) / ks) + 1 = (u32::MAX / ks) + 1 reciprocal: (u32::MAX / ks) + 1, rounding_bias: ks / 2, } } fn to_store(self, w: U8Weight) -> u8 { let n = self + w.rounding_bias; ((n as u64 * w.reciprocal as u64) >> 32) as u8 } } ``` ### 影响 (https://apas.tel/blog/optimizing-image-rs-blur#impact-1) 仅用倒数替换，就在整数累加器的基础上，对所有核大小获得了 ×3 的加速： fast_blur σ=50 ×2.8 更快（越低越好） 之前：整数加除法 之后：整数加倒数 这证实了除法是去掉浮点后内循环中的主要成本，而用倒数乘法替代是一个巨大的胜利。 ## 最终结果 (https://apas.tel/blog/optimizing-image-rs-blur#where-this-lands) 总体而言，整数累加器和倒数乘法的组合在 `u8` 图像上实现了 **5.9 倍的加速**，从约 52ms 降至约 8ms。这意味着每秒可以处理 120 帧而不是 19 帧，足够用于视频流或游戏等实时应用。 这两项优化都归结于一个观察：并非所有指令都平等。浮点指令的代价比整数指令高好几个数量级。`div` 的代价比 `mul` 加移位大约高一个数量级。可以确定的是，CPU 流水线非常迷人，我还没有探索完。 完整的实现已被合并 (https://github.com/image-rs/image/pull/2846)，并将在 image-rs 的下一个版本中可用。 感谢 @awxkee (https://github.com/awxkee) 提供算术专业知识，感谢 @197g (https://github.com/197g) 和 @RunDevelopment (https://github.com/RunDevelopment) 提供设计反馈。 1. Kovesi, P.: **Fast Almost-Gaussian Filtering.** 澳大利亚模式识别学会会议：DICTA 2010. 2010 年 12 月. 悉尼. PDF (http://www.peterkovesi.com/papers/FastGaussianSmoothing.pdf) ↩ (https://apas.tel/blog/optimizing-image-rs-blur#user-content-fnref-1) 2. Fog, A.: **指令表：Intel、AMD 和 VIA CPU 的指令延迟、吞吐量和微操作分解列表.** PDF (https://www.agner.org/optimize/instruction_tables.pdf) ↩ (https://apas.tel/blog/optimizing-image-rs-blur#user-content-fnref-2) 3. 基于 Granlund, T., Montgomery, P.L.: **通过乘法实现不变整数除法.** ACM SIGPLAN Notices, 29(6), 1994. PDF (https://gmplib.org/~tege/divcnst-pldi94.pdf) ↩ (https://apas.tel/blog/optimizing-image-rs-blur#user-content-fnref-3)