5x5浮点中值滤波优化实战:从std::sort到无分支排序网络
2026/9/16 4:39:29 网站建设 项目流程

前阵子接手一个激光雷达深度图去噪项目,数据是单精度浮点,要求做5x5中值滤波。一开始图省事,直接把每个窗口的25个浮点数丢进std::sort,排完取第13个,逻辑倒是清爽,上板一测直接傻眼——单帧耗时根本压不住实时要求。后来把“5x5浮点中值滤波”当成一个正经的算法优化题目来啃,从比较次数、指令选择到内存访问模式一点点抠,最终把单窗口耗时压到了原来的三分之一以下。这篇文章就是整个优化过程的完整记录,包括方案取舍、浮点数据特有的坑、实测数据,以及几个值得记下来的工程细节。如果你正在做嵌入式图像处理、信号滤波,或者需要在实时系统里处理浮点窗口数据,这篇应该能帮上忙。

1. 为什么5x5浮点中值滤波值得专门优化

先说说中值滤波本身。它的原理非常朴素:取当前像素周围一个n×n窗口,把所有像素值排序,取中间那个值作为输出。相比均值滤波,中值滤波最大的优势是能在抑制离群点和脉冲噪声的同时保留边缘信息——均值会把边缘抹糊,中值不会。这也是它在激光雷达深度图、红外图像、医学影像里被广泛使用的原因。

但5x5窗口和常用的3x3窗口相比,复杂度完全不是一个量级。3x3窗口是9个数找中值,5x5窗口是25个数找中值。很多人觉得“不就多几个数嘛”,但排序的比较次数是随数据规模平方增长的,25个数的全排序比较次数是9个数的七八倍。更麻烦的是,如果只用中值结果而不需要完整排序列表,那全排序本身就有大量浪费。

另一个关键点是数据格式。如果是8bit灰度图,可以做256个桶的直方图,滑动窗口时增删桶,效率极高。但浮点数据没有这种福利——float的取值空间太大,直方图方案直接失效,只能走“比较”这条路。而浮点比较在嵌入式平台上并不便宜,排序网络、快速选择这类纯比较方案就成了主战场。

做这个优化之前,我先明确了两个工程目标:

  • 每个输出像素的计算时间必须是确定性的,不能出现最坏情况退化。
  • 代码要能方便地移植到不同嵌入式平台,至少不能在编译器开优化之后反而变慢。

确定这三点之后,整个优化方向就很清晰:把“一堆数排序取中间”这个问题,改写成“固定次数、无分支、可并行”的指令序列。

2. 25个数找中值的复杂度账,先算清楚再动手

拿到一个优化题目,我习惯先把复杂度账算明白。特别是这种窗口大小固定的算法,比较次数直接决定了性能上限。

2.1 全排序为什么是个坑:300次比较的浪费在哪里

最暴力的方案,就是把25个float全部排序。用一个基于比较的排序算法,至少需要O(n log n)次比较,实际工程里更多。以25个元素为例,冒泡排序需要300次比较,快速排序平均也需要100多次比较,std::sort作为内省排序,在25这个规模上表现也不理想——除了比较本身,还会引入递归调用、栈操作、迭代器移动这些额外开销。

但问题是:我们真的需要完整的排序序列吗?中值滤波只需要第13个元素,也就是有序序列正中间的那个。为了这一个数,把其他24个数的顺序也都排出来,这是最大的浪费。就像你只需要知道全班第13名的成绩,却把全班的成绩单从头到尾排了一遍,这中间大部分排序操作对最终结果毫无贡献。

除此之外,全排序还有个问题——分支多。基于比较的排序算法内部充斥着if (a > b)这类条件跳转。在嵌入式CPU上,分支预测失败的代价是很高的,流水线会被冲刷掉十几个周期。25个元素的排序可能有几十次分支跳转,每次都可能预测失败,实际耗时远比“比较次数×单次比较耗时”要难看。

2.2 快速选择:平均很快,但硬实时系统不敢用

那不全排序,用快速选择(quickselect)行不行?思路是找第13小的元素,而不是全排序。算法会在平均情况下把比较次数压到100次左右,比全排序快不少。实测下来,在桌面平台上快速选择确实能跑出不错的成绩。

但这里有个核心问题:快速选择的分区操作基于枢轴值,最坏情况会退化成O(n²)。虽然25个元素的窗口内,退化的概率不高,但对一个需要硬实时保证的系统来说,最坏情况耗时是多少必须明确。我们做嵌入式滤波时,系统每帧的处理时间是有预算的,不能接受“这次运气不好多跑了三倍时间”这种事。

另一个问题是递归。标准实现是递归的,每次递归都要保存上下文、传递参数。在资源受限的MCU上,哪怕只是几层递归,栈和调用开销也不容忽视。当然可以手写迭代版,但代码复杂度就上去了。

所以快速选择适合对平均性能敏感、不要求最坏情况保证的场景。如果是硬实时或者高确定性要求的场景,有更好的选择。

2.3 排序网络:固定比较序列、无分支、天然可展开

排序网络是一种“旁路分支”的思路:我们不写if (a > b),而是写一个“比较-交换器”,保证左边的值永远小于等于右边。整个排序过程变成一串固定的比较-交换操作序列,不再依赖数据内容跳转。

这个方案有三个突出优点:

  • 比较次数固定,执行时间完全确定,可以做WCET(最坏执行时间)分析。
  • 代码里没有分支,不存在分支预测失败。
  • 整个网络可以直接展开成一条直线指令流,还可以用流水线方式执行——前一次的min/max结果可以直接作为下一次的输入。

排序网络的关键在于构建比较-交换序列。对25个输入,可以搜到一个深度合理的排序网络,也可以程序化生成。比较次数虽然比理论最优值略多,但远低于冒泡排序,而且无分支带来的收益在嵌入式平台上尤为明显。

3. 排序网络优化落地:把比较-交换变成两条浮点指令

排序网络的核心元件是“比较-交换器”。在交叉点上比较两个数,如果顺序不对,就交换位置。传统写法是:

if (a > b) { tmp = a; a = b; b = tmp; }

这段代码在嵌入式平台上有两个问题:一是有分支,二是要搬数据。但如果我们真的要优化的对象是float,就有更漂亮的写法。

3.1 比较-交换器就是min/max两条指令

对于浮点数来说,比较-交换器的本质是:把较小的值放到低位置,把较大的值放到高位置。这不就是fminfmax吗?ARM平台上对应fminnmfmaxnm这两条指令,x86平台上也有对应的vminpsvmaxps。用它们实现比较-交换,既没有分支,也不需要中间变量,还不需要条件跳转。

static inline void cmpswap(float *a, float *b) { float lo = fminf(*a, *b); float hi = fmaxf(*a, *b); *a = lo; *b = hi; }

fminf/fmaxf替换之后,编译器在开启优化时能把这段代码映射到硬件的min/max指令上。实测在ARM Cortex-M7上,整个比较-交换过程只用两条浮点指令就完成了,前面那段带分支的版本要花费包括比较、跳转、搬运在内的五六条指令,而且还有分支预测失败的风险。

这里有个细节值得注意:GCC在不开-ffast-math时,对fminffmaxf的代码生成有时会比较保守。如果发现编译器没把它映射到硬件指令,可以改用平台内置函数,比如ARM的vminnmq_f32、x86的_mm_min_ps,或者GCC的__builtin_fminf。我在项目里就是直接用__builtin_fminf/__builtin_fmaxf,保证代码生成可控。

3.2 用Batcher奇偶归并网络生成5x5的排序结构

那么问题来了:这个比较-交换序列应该怎么组织?手写25个输入的排序网络不太现实,容易出错。常规做法是用Batcher奇偶归并排序网络程序化生成。

Batcher奇偶归并的思路是:先把序列拆成两部分,各自排序,然后把两部分按规则交叉比较归并。对25这个规模,生成器输出的比较序列很规律,可以在编译期用模板或宏展开,也可以在初始化阶段生成好比较对数组并复用。

核心代码如下:

// 生成Batcher奇偶归并网络的比较操作序列(示意) void batcher_oddeven_merge(int lo, int n, int r) { int step = r * 2; if (step < n) { batcher_oddeven_merge(lo, n, step); batcher_oddeven_merge(lo + r, n, step); for (int i = lo + r; i + r < lo + n; i += step) { compare_and_swap(i, i + r); } } else { compare_and_swap(lo, lo + r); } }

当然,Batcher奇偶归并排序网络不是比较次数最少的排序网络,但它的优势是生成规则简单、代码短,嵌入式场景下“稍微多几次比较但结构规整”比“拼比较次数极限但代码复杂”更划算。而且在25这个规模上,Batcher网络产生的比较次数已经是三位数以内,配合min/max指令,性能完全够用。

还有一点:在实际工程里,我觉得值得把生成的比较序列导出来硬编码,而不是每次运行时调用生成器。这样既能避免生成器本身的指令开销,又能让编译器把整个网络展开成直线代码,充分发挥无分支优势。

3.3 不需要全排序:只求第13个元素的部分网络

如果继续深挖,排序网络其实还可以进一步做“部分网络”——只跟踪第13个元素的位置,不关心其他元素的最终顺序。这就是选择网络(selection network)的思路。

选择网络在结构上仍然是固定的比较-交换序列,但它只保证第k个位置的元素是正确的。相比全排序网络,比较次数可以进一步降低。对25个元素求中值,理论上可以找到一个比较次数明显少于全排序网络的选择网络。我在工程里没有手工推导这个网络,而是用脚本搜索了一圈,找到一组适合25输入的序列,然后硬编码进代码里。

这部分思路值得展开说:你完全可以写出一个小脚本,基于“只保留第k个位置结果”的约束去剪枝生成比较序列。剪枝的原理很简单——全排序网络里有些比较操作只影响第k个元素之外的位置,如果目标只是第k个元素,这些比较可能会被部分省略。生成出来的序列仍然是无分支的,而且比较次数比完整排序网络更少。

3.4 工程上的近似方案:分层中值,把代价再压一个量级

如果你对滤波结果的精度要求不是“严格中值”,而是“基本能滤掉离群点”,那还有一个更快的近似方案:分层中值。思路是把5x5窗口拆成5行,每行5个元素先排一次序,选出每行的中值,得到5个代表值;再对这5个代表值排序,取其中值作为最终输出。

这样做比较次数会大幅下降。实测很多场景下,这个近似中值的中值(Median of Medians)在深度图去噪里效果和精确中值几乎看不出差别,但速度能再快一倍左右。代价是它不再是数学意义上的中值,对某些数据分布可能会引入轻微偏差。

我通常的建议是:项目初期先用精确中值把流程跑通,性能不够时再评估能否切换到这个近似方案。毕竟“能实时运行的近似结果”远好过“卡成PPT的精确结果”。

4. 浮点数据特有的坑:NaN、负零与按位比较技巧

用排序网络处理浮点数据,我踩过的坑比预想的多。下面这几个问题,任何一个不注意都会导致结果错误,而且很难排查。

4.1 NaN会毁掉整个排序网络

IEEE 754标准里,NaN和任何数比较都是false。也就是说(NaN > a)为false,(NaN < a)也为false,(NaN == NaN)还是false。如果窗口里混入一个NaN,排序网络里所有的比较-交换器都会“不知道该拿它怎么办”,最终结果可能是任意值——既可能是NaN,也可能是某个完全不应该出现的数,取决于比较器的实现。

激光雷达深度图里出现NaN并不罕见:无效测距点、信号丢失、强反射区域,都可能产出NaN。处理办法有两类:

  • 在滤波前扫描窗口,如果发现NaN,统一替换成一个固定的极小值(比如-FLT_MAX),把NaN当成“最不可能被选中”的离群点处理。
  • 如果业务上要求NaN传播(即输出也应该是NaN),那就单独走一条分支,不再进入排序网络。

我选了第一种,因为我们的后处理管线里NaN本来就要被过滤,输出NaN只会让后续模块多写一套处理逻辑。替换成极小值之后,中值滤波天然会把它们排在后面,不会污染输出。

4.2 IEEE 754位模式反转:把浮点比较变成整数比较

有时候我们希望避免浮点比较指令,或者想把浮点数据塞进整数排序网络路径。经典的做法是:在排除NaN之后,可以通过一个位变换把浮点数映射为“顺序一致”的32位无符号整数,之后就能用整数比较器。

核心技巧如下:

  • 正浮点数:符号位是0,直接翻转符号位(变成1),使所有正数的映射结果大于所有负数的映射结果。
  • 负浮点数:符号位是1,对所有位取反,使得绝对值较大的负数映射为较小的整数。

代码实现:

static inline uint32_t f32_sort_key(float f) { uint32_t u; memcpy(&u, &f, sizeof(u)); // 如果符号位为负,取反全部位;否则只翻转符号位 uint32_t mask = (uint32_t)(-(int32_t)(u >> 31)) | 0x80000000u; return u ^ mask; }

这个技巧的优势在于:整数比较没有NaN干扰(前提是已经清掉NaN)、没有浮点异常标志位被设置的问题、在某些只擅长整数运算的DSP上速度更快。它的代价是多了一次位变换的开销。在ARM Cortex-M系列上,由于存在fminnm/fmaxnm这类硬件指令,位变换方案未必更优;但如果你碰到一个没有原生浮点比较的廉价MCU,这个方案就是救命稻草。

4.3 正零和负零,以及大量重复数据

还有一个容易忽略的细节:+0.0f-0.0f在IEEE 754里是不同的位模式,但比较时它们相等。如果窗口里同时出现正零和负零,排序网络对这种“相等但位模式不同”的数据处理没有问题,输出中值可能是两者之一,这对后续业务没有影响。

但如果你加上面的位变换再去比较,正零和负零会被映射到不同的整数位置:正零映射到0x80000000,负零映射到0x7FFFFFFF。如果用整数排序网络选出中值位模式,返回时直接memcpy回float,输出可能是-0.0f。这个差异在大多数场景下无害,但在一些严格校验输出的测试用例里会翻车。我后来在实现里做了个统一处理:位变换之前先把-0.0f归一化为+0.0f

重复数据的问题同样值得说。中值滤波的窗口里经常出现大量相同值,比如深度图里平坦区域。排序网络对重复值的处理是稳定的吗?答案是:只要比较-交换器能做到“相等时保持原顺序”,排序网络的结果就是确定的。虽然中值不受小扰动影响,但为了保证多帧之间的输出连贯性,我建议在比较-交换器里对严格大于才交换,等于时不交换。

5. 实测数据、平台差异,以及几个值得记录的优化细节

理论说了一堆,最后落实到数字上。这部分数据来自我手头一块Cortex-M7,主频400MHz,单精度浮点,开了-O2优化。每个数字都是工程上的量级参考,不同编译器版本、不同优化选项下会有差异,但趋势是稳定的。

5.1 几种方案在单窗口(25个float)上的实测对比

方案单窗口近似耗时比较次数(约)是否有分支最坏情况
std::sort全排序1.1 us120~300有大量分支固定但开销最大
手写快速选择0.4 us90~110有分支但较少O(n²)退化风险
Batcher排序网络+min/max0.33 us固定150~180无分支完全固定
选择网络求第13个元素0.25 us固定110~130无分支完全固定
近似分层中值0.15 us固定75无分支完全固定

从上到下,性能逐步提升,代价是代码复杂度和“数学严格性”逐步下降。我在项目里最终选了“选择网络求第13个元素”这个方案,因为它在精确中值的约束下做到了最快,而且执行时间固定,非常适合实时系统。

5.2 滑动窗口访存优化:别傻傻重复读25个点

性能优化的另一个大头是内存访问。处理整幅图时,如果每个输出像素都重新从原始数据里读25个float,相邻窗口之间会有大量重复访问。5x5窗口向右滑动一个像素时,其实只有最右边一列的5个点是新的,左边有5个点被丢弃,中间20个点上一轮已经读过了。

工程上的解法是行缓冲(line buffer)。维护最近5行的数据,滤波时每次只从原始数据里读入5个新浮点数,而不是25个。这样访存带宽需求直接降为原来的五分之一。对缓存很小的嵌入式平台,这个优化非常关键。

边缘处理也要提前想清楚。5x5窗口无法覆盖图像的四条边,常见的做法是复制边缘像素(clamp)或者直接跳过边缘输出。对浮点数据,我强烈建议不要用0填充,因为0在深度图里往往代表“无效距离”,会强行拉低窗口的数值分布,让中值偏移到错误的位置。我用的是clamp策略:越界坐标就用最近的边界像素值填充,虽然多写几行分支,但滤波质量有保障。

5.3 多窗口向量化:一次同时处理4个像素

如果目标平台有SIMD指令(ARM的NEON、x86的SSE/AVX),还有一个更激进的优化思路:同时处理多个窗口。

5x5窗口的结构对所有像素都是一样的,排序网络里的每一层比较-交换都是“位置上固定的两个数互换”。与其一次处理一个像素,不如把4个不同窗口的“同一位置”的数据打包成float32x4_t向量,然后对向量执行min/max操作。

这意味着排序网络代码几乎不用改,只需要把标量fminf/fmaxf换成NEON的vminnmq_f32/vmaxnmq_f32,每个比较-交换器一次就能同时完成4个窗口的计算。吞吐量提升是接近4倍的(受限于寄存器压力和内存带宽,实际能到2.5~3.5倍)。

但要注意一个问题:SIMD版本的排序网络需要25个向量寄存器来保存中间状态。ARMv7 NEON有32个64位寄存器(或16个128位寄存器),如果只需要存float32x4_t,25个向量用掉25个寄存器,加上临时变量,寄存器分配会非常紧张。编译器可能会溢出,反而变慢。我当时的解决办法是调整网络结构,尽量让寄存器的生存周期错开,或者干脆一次处理2个窗口而不是4个,让寄存器压力降下来。

5.4 内存对齐、编译选项和inline带来的差异

最后几个细节,都是实际跑代码才会遇到的。

内存对齐:如果要对浮点数组做NEON向量加载,数组首地址必须16字节对齐。我用alignas(16)来声明行缓冲数组,避免在运行时动态对齐导致的开销。

编译器选项:不开任何fast-math,因为中值滤波本身不依赖浮点数学重排,但我也没有开-ffast-math,主要担心影响代码里其他部分的浮点计算。对于这个文件,我会确认编译器把__builtin_fminf/__builtin_fmaxf正确映射到硬件指令。

函数内联:排序网络被展开后会产生大量的比较-交换调用,如果比较器不是inline,函数调用开销会吃掉所有优化红利。声明成static inline是最基本的操作。我还顺手加了__attribute__((always_inline)),确保网络展开成直线代码。

指针别名:如果滤波是in-place操作(输入输出同一块内存),记得给函数参数加上restrict关键字。否则编译器会保守地假设读写可能冲突,不敢把两次加载重排,性能会损失不少。这个问题藏得很深,一度让我误以为是排序网络本身的问题。

最后补一句

中值滤波这类窗口算法,性能瓶颈从来不在内存带宽,而在比较操作本身。每个输出像素固定读25个float,计算强度远高于访存强度。所以优化的核心就一句话:把比较次数压下来,把每次比较的代价压下来。排序网络负责前者,min/max指令和SIMD负责后者。这套思路不止适用于中值滤波,凡是滑动窗口内求分位数、做形态学滤波的场景,都可以直接套用。我做这个项目时最大的体会是:不要被“25个数排序取中间”这个朴素描述骗了,算法方案的差别在这种小规模数据上,远比想象中大得多。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询