因有關需求，最近要弄個5*5的浮點版本的中值濾波， 而且希望效率能做到極致。感覺這個事情不是隨手搞定嘛，因為5*5的字節版本的中值已經有模版了，改為浮點的大差不差的，沒有啥問題。
確實如此，直接把早期的自己版本代碼弄過來，改下相關的數據類型和函數，功能一下子就出來了，而且效果也不俗。
在【算法隨記三】小半徑中值模糊的急速實現（16MB圖7.5ms實現） + Photoshop中蒙塵和劃痕算法解讀 中，曾經提到對於5*5的中值，理論上需要131次比較，但是我現在找不到這個131的數據的來源和依據了。
這幾天看Opencv的代碼，在其medianFilter.cl以及median_blur.simd.hpp裏的都有類似的比較，具體如下：
這裏的op算子就是比較算子，一共有113處。
在https://github.com/ravkum/medianFilter/blob/main/src/medianFilter.cl中，也有類似的算子，如下所示：
這裏只有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)
因此，這樣處理實際上是個近似的結果。
以下是一些耗時的比較結果：
可見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個浮點數。
