因有关需求，最近要弄个5*5的浮点版本的中值滤波， 而且希望效率能做到极致。感觉这个事情不是随手搞定嘛，因为5*5的字节版本的中值已经有模版了，改为浮点的大差不差的，没有啥问题。
确实如此，直接把早期的自己版本代码弄过来，改下相关的数据类型和函数，功能一下子就出来了，而且效果也不俗。
在【算法随记三】小半径中值模糊的急速实现（16MB图7.5ms实现） + Photoshop中蒙尘和划痕算法解读 中，曾经提到对于5*5的中值，理论上需要131次比较，但是我现在找不到这个131的数据的来源和依据了。
这几天看Opencv的代码，在其medianFilter.cl以及median_blur.simd.hpp里的都有类似的比较，具体如下：
image image
这里的op算子就是比较算子，一共有113处。
在https://github.com/ravkum/medianFilter/blob/main/src/medianFilter.cl中，也有类似的算子，如下所示：
image
这里只有99个算子。
所以优化的第一步就是用这99个算子的方案代替原先我收集的131个的方案，这个就会有25%左右的提速了。
而且测试这99算子的结果和原先的131个以及opencv的113个算子的结果是一样的。
而我们在百度搜索5*5中值滤波快速算法，可以得到以下的信息：
5x5中值滤波的快速算法旨在减少标准排序方法的高计算复杂度，通过优化排序或利用数据结构来加速中值查找。以下是一些常见的快速实现思路：
行-列排序优化：先对5x5窗口的每一行进行排序，确保每行有序；然后对每一列排序，使列内元素有序。经过此处理，窗口被分为四个象限，中心元素通常接近中值，从而减少比较次数。
在https://github.com/shrtique/sv_axis_median_5x5/blob/master/sources_1/new/median_5x5_top_module.sv中，我们可以看到如下的一段文字说明：
复制代码
//This is simplified module of median filtration with 5x5 window.
//It has 3 stages: sorting by lines -> sorting by columns -> sorting main diagonal*
//* main diagonal is one that has the smallest value of "biggest" column and..
//.. biggest value of "smallest" column
//
// EXAMPLE
//
// [ 20, 10, 1, 6, 8 ] [ 20, 10, 8, 6, 1 ] [ 110, 72, 50, 48, 36 ] [ * , * , *, *, 36 ]
// [ 32, 5, 40, 100, 60 ] [ 100, 60, 40, 32, 5 ] [ 100, 60, 40, 32, 7 ] [ * , * , *, 32, * ]
// [ 14, 3, 22, 26, 68 ] --> [ 68, 26, 22, 14, 3 ] --> [ 80 , 30, 22, 14, 5 ] --> [ * , * , 26, * , * ]
// [ 110, 12, 30, 7, 11 ] [ 110, 30, 12, 11, 7 ] [ 68 , 26, 12, 11, 3 ] [ * , 22, * , * , * ]
// [ 72, 48, 36, 50, 80 ] [ 80, 72, 50, 48, 36 ] [ 20 , 10, 8 , 6 , 1 ] [ 20, * , * , * , * ]
//
// median = 26
//
//IMORTANT NOTE:
//This algorithm could make a mistake in one discret (one position): so if we range all pixels in upwards direction..
//.. the median value will be on position 13 (from 1 up to 25), but sometimes algorithm takes neighbour value (e.g. 12).
//We decided that this mistake can't dramatically decrease efficiency of suppressing impulsive noise on image
复制代码
也是一样的意思： sorting by lines -> sorting by columns -> sorting main diagonal
那这样对于5*5的数据，就是5次水平方向排序，5次垂直方向排序，还有一次对角线排序，一共是11次排序，而且每次排序都是对5个数据极性排序，至于是升序还是降序其实不重要。
这个时候我们再看前面的几个位置的排序次数，反倒是那个99次排序很有意义，他恰好是11的整数倍，即每次排序需要进行9次比较。
根据冒泡排序的一般过程，如下所示：
复制代码
int len = 5;
int i = 0; int j; int t;
for (i = 0; i < len - 1; i++)
{
for (j = 0; j< len - i - 1; j++)
{
if (a[j]>a[j + 1])
{
t = a[j];
a[j] = a[j + 1];
a[j + 1] = t;
}
}
}
复制代码
假定五个数据分别为p0/p1/p2/p3/p4，则按照上面的循环展开后可得5个数据的排序需要以下的比较操作：
OP(p0, p1); OP(p1, p2); OP(p2, p3); OP(p3, p4);
OP(p0, p1); OP(p1, p2); OP(p2, p3);
OP(p0, p1); OP(p1, p2);
OP(p0, p1);
一共是10次，也不是上面说的9次。
不过观察opencv的那一部分代码的前面一部分，我们把下标为0到4的部分提取到一起，能得到如下的操作组合：
OP(p1, p2); OP(p0, p1); OP(p1, p2); OP(p3, p4);
OP(p0, p3); OP(p2, p3); OP(p1, p4);
OP(p1, p2); OP(p3, p4);
这里只有9个算子，实际拿数据测试，他们确实是足以完成5*5的一个排序的，所以opencv的数据集合github的代码组合在一起好像就能解释了99这个数据了。
不过我们观察ravkum的那个算子的排布，感觉又不是这个样子，整体好混乱的，先不管，我们就按照刚才了解的sorting by lines -> sorting by columns -> sorting main diagonal处理看是否能有进一步的处理空间。
首先sorting by lines这个没有啥好优化的，必须完整的进行，下一部， sorting by columns就有发挥的空间了，我们看上面的列举的数据：
// [ 20, 10, 1, 6, 8 ] [ 20, 10, 8, 6, 1 ] [ 110, 72, 50, 48, 36 ] [ * , * , *, *, 36 ]
// [ 32, 5, 40, 100, 60 ] [ 100, 60, 40, 32, 5 ] [ 100, 60, 40, 32, 7 ] [ * , * , *, 32, * ]
// [ 14, 3, 22, 26, 68 ] --> [ 68, 26, 22, 14, 3 ] --> [ 80 , 30, 22, 14, 5 ] --> [ * , * , 26, * , * ]
// [ 110, 12, 30, 7, 11 ] [ 110, 30, 12, 11, 7 ] [ 68 , 26, 12, 11, 3 ] [ * , 22, * , * , * ]
// [ 72, 48, 36, 50, 80 ] [ 80, 72, 50, 48, 36 ] [ 20 , 10, 8 , 6 , 1 ] [ 20, * , * , * , * ]
sorting by columns的目的是为了获取对角线的那五个数，可以考到第一列那个20，实际上是第一列的最小值，而第五列的36实际是该列的最大值，因此，这两列是完全不需要进行完整的排序的，只要找到的最小和最大值就可以，这样就节省了很大的一部分计算了。
第二列和第四列其实就是找到第二个最小值和第二个最大值，这个也不需要进行完整的排序，我们观察的我们的冒泡排序的4行算子，第一行
OP(p0, p1); OP(p1, p2); OP(p2, p3); OP(p3, p4);
处理完成后其实p4就是最大值或者最小值了（取决于op的操作），这个时候我们只要求p0/p1/p2/p3中的最大值或者最小值就可以了，而不用继续进行后续的排序了。基本上节省了一半的排序操作。
第三列本身就是求中值的过程，这里也可以用类似上面的操作，只不过这个是要进行冒泡的前两行操作：
OP(p0, p1); OP(p1, p2); OP(p2, p3); OP(p3, p4);
OP(p0, p1); OP(p1, p2); OP(p2, p3);
此时p4是最大值，p3是次最大值，因此只要比较p0/p1/p2中的最大值就可以获得中间值了，这个减少的工作量就比较有限了。
进行这样的优化，测试比直接使用99次比较算子，可以大概有个10%左右的提速。
但是，我们测试发现，在有些情况下，这个优化后的结果和原始的结果似乎不是完全一致，而用严格的冒泡排序对25个数据进行排序后的结果表明，原始的方式的结果才是正确的，后面仔细研究发现，其实几个地方的文档都有提及这个事情，百度的AI的结果有这样一句话：
中心元素通常接近中值，注意这里的接近。
而shrtique的注释了也写到了：
This algorithm could make a mistake in one discret (one position)
因此，这样处理实际上是个近似的结果。
以下是一些耗时的比较结果：
image
可见Opencv还是蛮快的，不过opencv有一个限制，就是浮点的数据，中值只支持3*3的或者5*5的。 而且这个函数线程数对耗时没有影响，而我这里可以开多线程，基本上一个线程就能提高一倍（不大于物理核）。
说一个小技巧，在这个加载一行5个数据的过程中，对于浮点数的SSE优化，普通的写法是：
__m128 P0 = _mm_loadu_ps(First + X);
__m128 P1 = _mm_loadu_ps(First + X + 1);
__m128 P2 = _mm_loadu_ps(First + X + 2);
__m128 P3 = _mm_loadu_ps(First + X + 3);
__m128 P4 = _mm_loadu_ps(First + X + 4);
那么考虑到数据的特殊性，其实可以用下面的操作比较过多的内存读写：
复制代码
// V0 V1 V2 V3 V4 V5 V6 V7
__m128 P0 = _mm_loadu_ps(First + X); // V0 V1 V2 V3
__m128 P4 = _mm_loadu_ps(First + X + 4); // V4 V5 V6 V7
__m128 P1 = _mm_castsi128_ps(_mm_alignr_epi8(_mm_castps_si128(P4), _mm_castps_si128(P0), 4)); // V1 V2 V3 V4
__m128 P2 = _mm_castsi128_ps(_mm_alignr_epi8(_mm_castps_si128(P4), _mm_castps_si128(P0), 8)); // V2 V3 V4 V5
__m128 P3 = _mm_castsi128_ps(_mm_alignr_epi8(_mm_castps_si128(P4), _mm_castps_si128(P0), 12)); // V3 V4 V5 V6
复制代码
这样只用加载2次内存，其他的数据用CPU指令组合而成，而上面的指令看上去很多，实际上那个cast并不产生任何的物理指令，他只是语法糖，所以就只有一个_mm_alignr_epi8指令，注意，这个应用不可以直接扩展到AVX，因为AVX不是把
_mm_alignr_epi8直接扩展到8个浮点数。
