常量时间中值滤波CTMF：原理、复现与实时图像去噪性能实测

1. 从一次图像预处理翻车说起：为什么中值滤波的耗时忽高忽低

几年前我在做一个工业质检的小项目，需要对流水线上采集的灰度图做实时去噪。图像分辨率不算大，大概 1280×1024，帧率要求 30fps。当时想都没想就上了经典的中值滤波，窗口 5×5，结果在测试机上跑得好好的，一上产线就出问题——帧率像坐过山车，有时候稳在 30fps，有时候直接掉到 12fps，而且掉帧的时机完全没规律。

排查了两天才反应过来：中值滤波的耗时和图像内容强相关。图像里平坦区域多的时候，排序很快；一旦遇到纹理复杂、噪声密集的区域，排序的数据分布变得"难缠"，比较次数飙升，单帧耗时直接翻倍。这就是典型的数据依赖型时间复杂度——算法跑多快，取决于你喂给它什么数据。

这个问题在实时系统里是致命的。你没法跟产线说"今天产品表面花纹复杂，所以检测速度慢一点"。于是我开始找有没有一种中值滤波，不管图像内容长什么样，耗时都是固定的。后来就翻到了 CTMF 这篇论文——Constant Time Median Filtering，常量时间中值滤波。它的核心承诺很直接：滤波耗时只和图像尺寸、窗口大小有关，和图像内容完全无关。

这篇博文就把我对 CTMF 的研读和复现过程完整拆开。我会讲清楚它到底怎么做到"常量时间"、和传统中值滤波的本质区别在哪、复现时哪些细节容易踩坑、以及在实际工程里它到底值不值得用。如果你也在做实时图像处理、嵌入式视觉或者对算法复杂度优化感兴趣，这篇内容应该能帮你省下不少自己啃论文的时间。

2. 传统中值滤波的复杂度到底卡在哪

2.1 朴素实现的 O(k² log k) 是怎么来的

先把这个基础问题说透，不然后面理解 CTMF 的巧思就没有参照系。

中值滤波的操作本身很朴素：对图像上每个像素，取以它为中心的 k×k 窗口，把这 k² 个像素值排序，取中间那个值作为输出。窗口 3×3 就是 9 个数排序取第 5 个，5×5 就是 25 个数排序取第 13 个。

如果每个像素都独立地重新排序，单像素的复杂度是 O(k² log k)，整幅图 W×H 就是 O(W·H·k² log k)。这个复杂度是固定的，和图像内容无关。但问题在于，没人会真的每个像素都从头排序——那样太慢了，工程上根本不可接受。

2.2 滑动窗口优化：把复杂度降下来，却引入了数据依赖

实际实现里，大家都会利用窗口的滑动特性做优化。当窗口从左边滑到右边时，只有最左一列移出、最右一列移入，中间大部分数据是复用的。基于这个观察，出现了很多优化结构：

• 直方图法：维护窗口内像素值的直方图，滑动时只更新移入移出的像素。对于 8 位图像，值域只有 0~255，直方图数组大小固定 256。每次滑动更新 2k 个计数，然后从直方图里找中位数。
• 排序树/平衡树：维护一个有序结构，插入删除都是 O(log k²)，取中位数 O(1) 或 O(log k²)。
• 堆结构：用两个堆维护中位数，插入删除 O(log k²)。

这些方法把单像素的均摊复杂度降到了 O(k) 甚至 O(log k)，看起来很美。但注意，"找中位数"这一步的耗时是数据依赖的。

拿直方图法举例：更新完直方图后，你要从 0 开始累加计数，直到累加值超过窗口元素总数的一半，那个位置就是中位数。如果中位数靠近值域中间，你要累加一百多次；如果图像整体偏暗，中位数靠近 0，累加几次就找到了。图像越亮或越暗，找中位数越快；图像灰度分布越均匀，越慢。这就是我当年产线上掉帧的根源。

这里有个容易被忽略的点：很多人以为直方图法已经是 O(1) 了，因为值域固定 256。但严格说它是 O(值域)，值域固定时算 O(1)，可这个"1"是个不小的常数，而且实际执行的循环次数随数据变化。在硬实时场景里，常数大不是问题，常数不稳定才是问题。

2.3 数据依赖为什么在实时系统里是原罪

在普通应用里，数据依赖导致的耗时波动可能无所谓，平均下来能接受就行。但在实时系统里，判断标准不是平均值，而是最坏情况执行时间（WCET，Worst-Case Execution Time）。

假设你的系统每帧预算 33ms（30fps），传统中值滤波平均耗时 20ms，但最坏情况 45ms。那你在设计时就必须按 45ms 来留余量，否则遇到最坏情况就崩。这意味着你白白浪费了 25ms 的算力预算，或者被迫降分辨率、降帧率。

CTMF 的价值就在这里：它把最坏情况拉平到和平均情况一样，让你可以精确地做时间预算。你不需要为"万一图像很复杂"预留额外余量，这对嵌入式设备和硬实时系统来说是实打实的收益。

3. CTMF 的核心思路：把"找中位数"变成"数数"

3.1 一个反直觉的切入点：中位数不需要排序

理解 CTMF 的第一个关键，是意识到求中位数根本不需要完整排序。

中位数的定义是：在 n 个元素中，有至少 n/2 个元素小于等于它，且至少 n/2 个元素大于等于它。换句话说，你只需要知道"有多少个数比某个候选值小"，就能定位中位数，不需要知道这些数的具体顺序。

这个观察把问题从"排序问题"转化成了"计数问题"。而计数问题，恰恰是可以做到和数据分布无关的——只要你换一种计数方式。

3.2 列直方图：CTMF 的数据结构基石

CTMF 论文里最核心的数据结构叫列直方图（column histogram）。这是它区别于传统直方图法的关键。

传统直方图法维护的是整个窗口的直方图，一个 256 大小的数组，滑动时更新移入移出的 2k 个像素。而 CTMF 维护的是每一列的直方图。

具体来说，对于一个 k×k 的窗口，它有 k 列。CTMF 为每一列单独维护一个直方图，记录这一列中 k 个像素值的分布。所以总共维护 k 个直方图，每个直方图大小等于值域（比如 256）。

当窗口向右滑动一列时：

• 最左边那一列的直方图被丢弃（或者说，它的角色被最右边新列取代）
• 新进入的那一列，需要重新统计它的直方图

等等，这里有个关键问题：如果每次滑动都要重新统计一整列的直方图，那复杂度不就上去了吗？别急，这正是 CTMF 精妙的地方，也是它和传统方法最大的不同。

3.3 为什么"每列独立统计"反而更快

传统直方图法滑动时，更新的是"移出的 k 个像素"和"移入的 k 个像素"，共 2k 次计数修改。CTMF 滑动时，看起来要重建一整列的直方图，也是 k 个像素，似乎没占到便宜。

但差别在于后续的查询操作。

传统方法更新完直方图后，要从头累加找中位数，这个累加次数是数据依赖的。CTMF 则不同，它维护的列直方图有一个特殊性质：每一列的直方图是"静态"的——在窗口停留期间，这一列的像素值不变，直方图也不变。只有整列移入移出时才更新。

这意味着 CTMF 可以把"找中位数"的操作，转化为对 k 个列直方图的并行累加。而且这个累加过程，可以通过预计算前缀和的方式，把每次查询的复杂度压到和数据无关。

论文里给出的做法是：对每一列的直方图，预先计算它的累积分布。这样查询"这一列有多少个像素小于等于 v"就是 O(1) 的。然后对 k 列做累加，找到中位数所在的列和值。整个过程是 k 次 O(1) 操作，和数据分布完全无关。

3.4 常量时间的真正含义：和数据无关，和窗口大小有关

这里必须澄清一个常见误解。CTMF 说的"常量时间"，指的是对每个像素的处理时间与图像内容无关，而不是说它和窗口大小 k 无关。

实际上，CTMF 的单像素复杂度是 O(k)，因为要对 k 列做累加。窗口越大，耗时越长。但关键是，这个 O(k) 是稳定的，不管图像是纯黑、纯白还是随机噪声，耗时都一样。

对比一下：

所以 CTMF 的定位很清晰：它用略高于最优均摊复杂度的代价，换来了最坏情况的确定性。在实时系统里，这笔买卖通常很划算。

4. 复现 CTMF 时最容易踩的四个坑

论文读起来逻辑很顺，但真正动手写代码时，魔鬼全在细节里。我复现的时候踩了不少坑，这里挑四个最有代表性的讲。

4.1 列直方图的更新时机搞错，导致结果偏移

第一个坑最隐蔽。我一开始以为，窗口每滑动一个像素，就要更新所有列的直方图。结果写出来的滤波结果整体偏移了一个像素，边缘还出现了奇怪的条纹。

问题出在对"列"的理解上。CTMF 的列直方图，是以窗口为单位维护的，不是以图像列为单位。当窗口向右滑动一个像素时，窗口覆盖的列整体右移一列，最左列移出、最右列移入。但窗口内部的列与列之间的相对关系没变。

正确的做法是：维护一个循环队列来存 k 个列直方图。滑动时，把最旧的列直方图重置，用新进入的那一列数据重新填充，然后把它放到队列尾部。这样队列里始终是当前窗口的 k 列，且顺序正确。

我当时的错误是每次滑动都重建所有列，不仅慢，而且因为重建顺序和窗口列顺序对不上，导致中位数取错。

4.2 累积分布预计算：别忘了处理边界

CTMF 查询中位数时，需要对每列的直方图做累积。我最初图省事，每次查询都现场累加，结果复杂度又退化成数据依赖了——因为累加次数取决于中位数落在值域的哪个位置。

正确做法是预计算累积分布。但这里有个边界问题：累积分布的最后一个元素应该等于该列的总像素数 k。如果你在更新直方图后忘记同步更新累积分布，或者累积分布的数组大小和值域对不上，就会出现"找不到中位数"的情况——累加到最后一个值还没超过总数的一半。

我的建议是，把直方图和累积分布封装成一个结构体，任何修改直方图的操作都强制同步更新累积分布。虽然多了一点开销，但能避免大量隐蔽 bug。

4.3 值域假设：8 位图像和 16 位图像的处理差异

论文里的例子默认是 8 位灰度图，值域 0~255。但实际项目里，我遇到过 16 位深度的医学图像，值域 0~65535。这时候如果直接套用 256 大小的直方图，就会溢出。

处理高 bit 深度图像时，有两个选择：

• 直接扩大直方图：用 65536 大小的数组。内存占用大，但逻辑简单。
• 量化降位：先把 16 位量化到 8 位再做滤波。会损失精度，但速度快。

我当时的项目对精度要求高，选了第一种。代价是每个列直方图占 64KB，k=5 时就是 320KB，对嵌入式设备来说内存压力不小。后来我改用了稀疏直方图——只记录实际出现的值，用哈希表存。这样在值域大但实际取值集中的场景下，内存和速度都能兼顾。

这里给个经验值：如果值域超过 4096，且图像实际灰度分布比较集中，优先考虑稀疏结构；如果分布很散，还是老老实实用稠密数组，哈希表的开销反而更大。

4.4 多线程下的数据竞争：列直方图不是只读的

最后一个坑和并行化有关。CTMF 的列直方图在窗口滑动时会被修改，如果你想把图像分块并行处理，块与块之间的边界列会被两个线程同时访问。

我当时的做法是让每个线程独立维护自己负责区域的列直方图，块边界处多算几列作为重叠区，最后合并时只取有效部分。这样虽然有一点重复计算，但避免了锁竞争。实测下来，4 线程能跑到接近 3.5 倍的加速比，比加锁方案快得多。

5. 实测对比：CTMF 到底比传统方法快多少

光讲原理不够，得拿数据说话。我在一台普通开发机上做了一组对比测试，图像尺寸 1920×1080，灰度图，分别测了 3×3、5×5、7×7 三种窗口。

5.1 测试环境与实现细节

测试用的是 C++ 实现，编译器开 O2 优化。对比对象有三个：

• 朴素排序法：每个像素独立排序，作为基准。
• 传统直方图法：维护窗口整体直方图，滑动更新。
• CTMF：列直方图 + 累积分布预计算。

测试图像选了三张，代表不同数据分布：

• 平坦图：大面积渐变，灰度分布集中。
• 纹理图：细节丰富，灰度分布均匀。
• 噪声图：随机噪声，灰度分布接近均匀但更散。

每张图跑 100 次取平均，同时记录单次最大耗时。

5.2 耗时数据：稳定性的价值一目了然

先看 5×5 窗口的结果：

这组数据把问题说得很清楚。传统直方图法在平坦图上最快，85ms，但一到噪声图就飙到 158ms，波动接近一倍。CTMF 虽然均值比传统方法在平坦图上慢一些（118ms vs 85ms），但三张图几乎一样快，最大耗时也只有 123ms。

在实时系统里，你要按最大耗时做预算。传统方法要按 165ms 留余量，CTMF 只要按 123ms。算下来 CTMF 反而能让你把帧率做更高，或者留出更多算力给后续处理。

5.3 窗口越大，CTMF 的优势越明显

再看不同窗口大小的对比（噪声图，最坏情况）：

规律很明显：窗口越大，传统直方图法的数据依赖越严重（因为要找的中位数位置更不确定），而 CTMF 的耗时增长接近线性，且始终稳定。到 9×9 时，CTMF 的最坏耗时已经不到传统方法的一半。

这也符合理论分析：CTMF 是 O(k)，传统直方图最坏情况接近 O(值域)，窗口越大，值域固定的劣势越明显。

5.4 什么场景该选 CTMF，什么场景不必

数据摆出来，选型逻辑就清晰了：

• 硬实时系统、嵌入式设备：闭眼选 CTMF。稳定性就是生命线，多花的那点平均耗时完全值得。
• 离线批处理、对时间不敏感：传统直方图法够用，实现简单，平均更快。
• 窗口特别大（7×7 以上）：CTMF 优势明显，即使非实时场景也建议用。
• 值域特别大（16 位以上）：要权衡内存，CTMF 的列直方图内存占用是 k 倍值域，需要评估。
• 图像内容高度可预测：比如固定场景的监控，传统方法可能更划算。

6. 把 CTMF 用进实际项目的几个工程建议

理论清楚了，数据也有了，最后聊聊落地。这部分是论文里不会写、但实际项目里一定会遇到的东西。

6.1 内存布局：列直方图怎么存才缓存友好

CTMF 要维护 k 个列直方图，每个直方图是一个数组。如果按"直方图1的所有元素、直方图2的所有元素……"这样顺序存，查询时对 k 个直方图做累加，会频繁跳内存，缓存命中率低。

更好的做法是转置存储：把所有列直方图的同一个值域位置放在一起。比如值 v=0 时，k 个列的计数连续存放；v=1 时又是连续 k 个。这样查询中位数时，从 v=0 开始逐值累加，每次访问都是连续内存，缓存友好度大幅提升。

我实测这个改动能带来 15%~20% 的性能提升，而且实现起来只是改一下索引计算方式，性价比很高。

6.2 边界处理：图像边缘的窗口怎么算

图像边缘的像素，窗口会超出图像范围。常见处理有三种：补零、复制边缘、镜像。CTMF 对这三种都支持，但要注意补零会改变直方图分布，可能影响中位数。

如果项目对边缘质量要求高，建议用镜像填充，它引入的伪影最少。如果只是内部区域做检测，边缘直接跳过最省事。我当时的产线项目就是只处理内部区域，边缘单独用简单方法处理，省了不少事。

6.3 和硬件加速的结合：SIMD 能不能用

CTMF 的累加操作天然适合 SIMD。k 个列的累积分布累加，可以打包成向量一次算。我用 SSE 试过，在 k=5 时能再快 30% 左右。不过要注意，SIMD 对数据对齐有要求，列直方图的内存布局要提前规划好。

如果目标平台有 GPU，CTMF 也能映射过去。每个像素的滤波独立，适合并行。但列直方图的维护有依赖，需要仔细设计 kernel 的划分方式。这块我还没深入做，有经验的读者可以补充。

6.4 一个容易被忽视的优化：提前终止

虽然 CTMF 号称常量时间，但在找中位数时，其实有一个可以提前终止的点：当累加计数已经超过窗口元素总数的一半时，就可以停了，不需要累加完所有值。

这个优化会让耗时重新变得和数据相关，但它只会让耗时变短，不会变长。也就是说，最坏情况仍然是常量时间，平均情况更快。在对最坏情况有要求、又想要平均性能的场景，这是个不错的折中。

我一般会保留这个优化，因为它不破坏 CTMF 的核心承诺——最坏情况可控。

7. 写在最后：常量时间思维的价值

复现完 CTMF 之后，我最大的收获其实不是这个算法本身，而是它背后的思维方式：当平均性能已经够好，但最坏情况不可控时，主动用一点平均性能去换确定性，往往是一笔划算的买卖。

这个思路在很多领域都适用。比如内存分配器设计、网络协议的超时处理、任务调度算法，都会遇到"平均快但最坏慢"和"平均稍慢但稳定"的取舍。CTMF 给了一个很漂亮的示范：通过换一种数据结构（列直方图替代整体直方图），把数据依赖从算法里彻底剥离出去。

如果你正在做实时图像处理，我建议至少把 CTMF 实现一遍，哪怕最后不用，这个过程也会让你对"复杂度"和"实时性"的理解深一个层次。代码我建议从 8 位灰度、3×3 窗口开始，跑通之后再扩展到更大窗口和更高位深。踩坑是必然的，但每踩一个，你对这个算法的理解就扎实一分。

后续如果要做视频序列的滤波，还可以考虑帧间复用——相邻帧的列直方图有大量重叠，理论上能进一步降低开销。这个方向我还在试，有结果再分享。