基于ARM架構的中值濾波算法優化*

牟明任，賈海鵬，張云泉，鄧明森，曲國遠，魏大洲，張廣婷

(1.貴州財經大學信息學院，貴州 貴陽 550025；2.中國科學院計算技術研究所計算機體系結構國家重點實驗室，北京 100190;3.中國航空無線電電子研究所，上海 200241)

1 引言

在圖像處理中，降低椒鹽噪聲的常用方法是中值濾波[1,2]。中值濾波是一種基于統計學的非線性濾波技術，圖像的噪聲值被中值濾波窗口內的中值所代替。中值濾波窗口遍歷整幅圖像，計算濾波窗口內所有值的中值作為新的像素值。中值濾波算法的中值計算公式如式(1)所示：

g(x,y)=median{f(x-i,y-i),i,j∈H×W}

(1)

其中,f(x,y)和g(x,y)分別是初始圖像的值和輸出圖像的替代值，H×W是濾波窗口的大小(通常H=W且為奇數，比如3×3，5×5，7×7…等)。

對于中值濾波的研究已經有很多，但是目前中值濾波算法仍存在一些問題。基于統計直方圖的中值濾波算法雖然在8u的灰度圖像的大濾波窗口下執行效率較高，但在小濾波窗口下并不適用，且只局限于在灰度圖像上使用。另一方面，針對非8u數據類型的大濾波窗口的中值濾波問題，目前還沒有很好的解決算法。

另外，在圖像處理問題上，SIMD(Single Instruction Multiple Data)指令是非常適用的，可以顯著提高圖像處理性能。然而之前的相關工作中并沒有充分利用這一并行方式。

本文針對不同數據類型以及不同規模濾波窗口下的自適應中值濾波算法,根據不同的算法提出了不同的優化策略。對于能夠并行的算法，本文使用ARM NEON指令集(ARM架構下的矢量寄存器可以以指令流水線的方式高效執行這些數據)進行SIMD實現，進一步提高圖像處理性能。

實驗結果表明,本文實現的HMPP算法庫中的中值濾波自適應算法帶來了較好的性能提升，相比于OpenCV的中值濾波算法，性能提升達到120%。

本文的主要貢獻如下：

(1)提出了中值濾波自適應算法，針對不同的濾波窗口大小及不同數據類型提出了相對應的算法，提高了中值濾波的效率。

(2)設計了一種針對ARMv8架構的優化策略，使用ARM NEON指令集進行優化，也可對其他算法在ARM處理器上的優化提供參考。

(3)實現了中值濾波算法，并與OpenCV進行了性能比較，且性能更好。

2 相關工作

現在已有較多關于中值濾波算法的研究工作[3-11]，基于分治策略的快速算法[5,6]將傳統中值濾波算法的時間復雜度從O(scale2)降低到O(scalelogscale)(其中scale表示圖像規模);對于8u數據類型的圖像來說，文獻[7]提出的算法基于可修改的直方圖——統計直方圖，有效地提升了中值濾波的效率。

在8u的圖像中，r，g和b的取值在0～255,那么就可以采用統計方式來計算濾波窗口內的中值，直方圖的數組長度為256，數組下標0～255的值為0～255。使用可滑動的濾波窗口，改變窗口內的像素值，滑動求取中值。一個像素的直方圖在窗口第一次運行時被設置，從一個位置滑動到下一個位置時，窗口內的值會改變，直到濾波窗口遍歷完整幅圖像。窗口的中值是通過使用更新的直方圖調整前一個窗口的中值來找到的。

對于處理灰度圖(8u數據類型)，當濾波窗口的半徑非常大時，使用直方圖算法是十分有效的，算法復雜度為O(scale)，且不會隨著濾波半徑的大小變化產生太大的性能波動。

之后研究人員又提出了一系列基于直方圖的改進算法。Perreault等[9]提出的改進算法，提高了中值濾波的效率，但是直方圖具有局限性，只能用于數據類型為8u的灰度圖像，對于具有更高位數的數據類型(16u,16s,32f,64f)的圖像來說，使用直方圖求取中值會產生巨大的內存開銷。以16u為例，一幅圖像中各個點的像素值在0～216,如果使用直方圖進行中值濾波,空間開銷極大，實現起來是不現實的，同時對于較小濾波窗口來說性能不佳。文獻[10]提出的算法中，將中值濾波與均值濾波相結合，能更好地抑制噪聲，保留圖像細節。引入依舊適用的統計直方圖來提高中值的搜索速度，充分利用了圖像的相關性。因此，改進算法的復雜度降低到了O(scale)。該算法的性能還取決于用于表示數據值的數據類型,在數據類型變化時，比較的規模以指數級增長。

3 自適應中值濾波算法

給定一個M×N大小的圖像，一個H×W(H

Figure 1 Finding the median value in the filtering window

由于中值濾波是一種非線性濾波，對于含有隨機噪聲的圖像，其數學分析較為復雜，對于均值濾波噪聲為零的正態分布圖像，中值濾波的噪聲方差[11]近似如式(2)所示：

(2)

(3)

比較式(2)和式(3)，中值濾波的性能依賴于濾波窗口的大小，同時由于圖像數據類型的不同，單一的算法并不能解決所有的問題，因此本文提出了自適應中值濾波算法，根據濾波窗口的大小和數據類型選擇最佳算法。

3.1 面向小濾波窗口的中值濾波通用算法

當濾波窗口較小時，即濾波窗口為3×3或5×5時，本文提出了一種單獨處理的算法。以3×3濾波窗口為例(H=3,W=3)，設濾波窗口內的像素值分別為P0、P1、P2、P3、P4、P5、P6、P7和P8。

本文使用歸并排序(Mergesort)的分組思想來優化算法。歸并排序算法首先將輸入序列分成若干個組，進行組內排序；然后對這些已排序的子序列再次分組，同組內的各有序序列歸并為一個有序序列，重復這個過程直到所有組合并為一個有序序列時算法結束。對于一個3×3的窗口，采用嚴格意義上的歸并算法得到中值需要11個時鐘周期。而此時所有數據均已排序，由于中值濾波只需找到中間值即可，因此不需要那么多次比較。首先，對每一行的3個像素進行排序，如圖2a所示，實現行方向的有序排序；然后將3組有序序列按其列有序進行排序，如圖2b所示。在最小值數據組中，1和2不可能是中值，因為窗口中至少有5個數據比它們大；在最大值數據組中，8和9不可能是中值，因為窗口中至少有5個數據比它們小；在中值數據組中，有4個以上的數據不能確定比它們大或者小，這樣3、4、5、6、7這5個數中的中間值就是最后窗口的中值。5個元素求中值需要6次比較，但是4、5、6已經有序，故節省了2次比較，可知最終通過16次比較就可以找到中值元素。如圖3c所示，取出有序序列組中的最小值中的最大值、最大值中的最小值和中間列的中值，3個數進行比較，就可得到中值。該算法的操作步驟見圖2，P0～P8分別代表3×3濾波窗口內的像素值，2個像素之間橫線下面的數字代表執行步驟。

Figure 2 Schematic diagram of median algorithm based on 3×3 filtering window

以圖2a中第①步的P1和P2為例，首先進行2個像素值的比較，當P1>P2時,P1和P2的像素值需要互換，否則繼續進行下一組像素值的比較。當濾波窗口大小為5×5時，算法步驟與3×3算法的一致。該算法可以應用于所有數據類型的3×3和5×5濾波窗口的中值濾波，算法復雜度為O(scale),同時可以通過ARM NEON指令集進行SIMD加速。

3.2 面向8u的中值濾波算法

當圖像為8u的灰度圖像時，min(W,H)≥7或者W≠H，本文會采用改進的直方圖算法進行中值濾波。

算法具體步驟如下：首先,構建初始直方圖，對于灰度圖像來說，像素的取值在0～255，設置一個數組，數組長度為256，從原始圖像數據(src)中獲取濾波窗口內的像素值，構建直方圖，獲取第一次濾波窗口內的像素值，對應的數組下標內的值加1，直到濾波窗口內的所有像素值全部放入數組內;然后,針對濾波窗口的大小設置閾值T=[(W×H)/2],依次累加數組內包含濾波窗口區域的值對應下標的個數sum，若sum≥T，就說明當前灰度層級的值為所需要的中值，同時濾波窗口向右滑動，依次采取上述方式進行濾波，遍歷整幅圖像。如圖3所示，虛線框與實線框內分別是前后2次濾波窗口內的數據。

Figure 3 Median filtering schematic diagram based on histogram

3.3 面向非8u的大濾波窗口的中值濾波算法

在數據類型的位數更多(16u,16s,32f,64f)的圖像中，直方圖根本無法使用，因此針對位數更多的圖像，要采取更加適合的算法，本文提出了快速中值算法(Fast Median Algorithm)，如算法1所示。

算法1快速中值算法

輸入：濾波窗口內所有的像素值。

輸出：所求得的中值。

步驟1獲取數列長度length=H×W,設有集合p(i),i∈{0,1,2,…,length-2,length-1}。

步驟2分別獲取i=0，(length-1)/2,length-1時的p(i)值，比較3個位置像數值的大小，最大值放在i=length-1，最小值放在i=(length-1)/2,i=0位置的像素值為哨兵，將i=1與i=(length-1)/2位置上的像素值互換。

步驟3通過與哨兵的像素值進行比較，確定最終i=m的位置，使得p(left)M(left∈{0,m-1},right∈{m+1,length-1})。

步驟4比較m與(length-1)/2的大小，如果m<(length-1)/2,更新數列為m～(length-1)；如果m>(length-1)/2,更新數列為1～m,繼續執行步驟1～步驟4。

4 中值濾波的優化策略

本文根據濾波窗口的大小采取最優的算法后，會使用ARM NEON指令進行SIMD優化。雖然編譯器在程序運行時會自動進行優化，但是根據研究發現[12],使用ARM NEON指令集進行手動優化，可以獲得更大的加速空間。在圖像處理中，使用SIMD可以發揮出更大的性能優勢，在ARM處理器中，寄存器被視為統一數據類型的元素的矢量,ARM的矢量寄存器擁有128位，向量寄存器分別可以設置為8，16，32和64寬度的位寬。以8u圖像為例，對于1個矢量寄存器1次可以存儲16個數據。

4.1 面向小濾波窗口的中值算法優化

對于3×3,5×5濾波窗口的中值濾波來說，本文已經通過3.1節的算法來降低時間復雜度，并且該算法可以使用ARM NEON指令進行優化。

單通道8u的灰度圖像的優化策略為：寄存器可以1次并行處理16個數據，意味著寄存器可以從原來1條指令比較2個數到現在1條 指令可以操作2個寄存器中的32個數據相比較一次性得出16個結果。首先使用vld1q_u8()將數據讀入矢量寄存器中。由于本文需要頻繁地比較數據大小，因此要使用vmin_u8()與vmax_u8()這2條指令，其中vmin_u8()是將2個數據中小的數據存入目標寄存器，vmax_u8()是將2個數據中大的數據存入目標寄存器，2條指令執行的操作都是對2組16個數據進行大小比較，并將結果存入到目標寄存器中。使用這2條指令按照3.1節提出的算法進行操作，最終執行1次可以得到16個中值，如圖4所示。

Figure 4 Register storage and operation for grayscale image

當8u的灰度圖像為多通道時，優化的難點在于，單通道是由r,g和b這3個像素中的一個決定的，有(r,r,r,…,r,r,r),(g,g,g,…,g,g,g),(b,b,b,…,b,b,b)3種形式，像素間沒有數據性依賴且連續，因此在濾波窗口滑動取值時比較容易讀入寄存器。

由于多通道灰度圖像由(r,g,b,…,r,g,b)或(r,g,b,a,…,r,g,b,a)組成,以第1個元素r為例，首先讀取第1個位置i，第2個位置為i+3。如果不采用矢量寄存器進行并行優化，該取數問題非常好解決，只需要控制數組下標即可，但是矢量寄存器讀入的數據是連續的，因此本文需要使數據按照需求讀入。

因此，本文提出以下優化策略：以8u圖像3通道為例，首先使用memcpy()函數將48個384位的r,g,b數據讀入內存中，r,g,b數量各為16個，但一個矢量寄存器的位數為128，那么本文需要3個這樣的矢量寄存器;使用vld3q_u8()代替uint8x16_tvld1q_u8()，將3組128位的數據依次讀入寄存器，這樣就會得到連續的r,g和b,并將其分別存儲在3個矢量寄存器中。

如圖5所示，本文通過多次使用vmin_u8()和vmax_u8()函數進行中值運算。求取中值之后采用vstd3q_u8()函數依次將數據存入目的圖像位置。該優化策略采用SIMD進行加速，理論上可以將性能提升16倍。同時，由于本文采取的算法將時間復雜度由原來的O(scale2)降低為O(scale),能很大幅度地提升小濾波窗口中值算法性能。本文在處理大規模數據時，性能十分穩定，在工程中具有較高的實用價值。

Figure 5 Multichannel register storage

4.2 基于8u直方圖改進的中值濾波優化

對于8u灰度圖min(W,H)>5且W≠H時，本文使用基于直方圖的中值濾波算法進行處理。在使用直方圖進行中值運算時，因為濾波窗口是以滑動的方式進行濾波，會有大量的數據被多次使用。本文在構建直方圖的數組時，會剔除垃圾數據，將新數據寫入數組中，但這樣會增加數據的復用性，以3×3濾波窗口為例，2個相鄰數據進行圖像濾波時，當前一個位置濾波后，濾波窗口滑動到后一個位置，會產生6個重復數據，如圖6所示。如果對大規模數據進行處理，如果數據的復用性不強，會造成資源的浪費，同時也會降低中值濾波算法的性能。根據本文分析，當濾波窗口從前一個位置滑動到后一個位置時，只需對直方圖進行一列加操作和一列減操作。如圖6所示，虛線窗口為前一個位置像素構成的濾波窗口，實線為后一個位置像素構成的濾波窗口。該算法的執行步驟如算法2所示。該算法可以增強數據的復用性，縮短數據讀取的時間，有效地減少運算量。該算法的優勢在于，當濾波窗口尺寸設置得很大時，該算法的濾波時間不會有很大的增加。

Figure 6 Median filtering optimization based on histogram

算法2基于直方圖的中值濾波算法

輸入:大小為M×N的圖像X，濾波窗口radio。

輸出:與圖像X大小相同的圖像Y。

//初始化直方圖H

//i,j分別為圖像X數組的橫縱坐標索引

//k表示濾波半徑內像素的索引

1：fori=1 toMdo

2:forj=1 toNdo

3:fork=-radiotoradiodo

4: RemoveXi+k,j-radio-1fromH;/*移除左邊一列*/

5: AddXi+k,j+radiotoH;//增加右邊一列

6:endfor

7:Yi,j←(median(H);/*median()函數對處理好的直方圖進行尋找中值的操作*/

8:endfor

9:endfor

當更新直方圖時，減掉當前濾波窗口的左側一列，加上當前窗口的右側一列，解決了數據重復使用的問題，開始對數據比較進行優化。更新完的新直方圖中，sum的值與T的關系有3種:當sum=T時,其中值維持不變；當sum>T時，說明中值在當前直方圖灰度層級的左邊;當sum

算法3基于改進直方圖的中值濾波算法

輸入:大小為M×N的圖像X，濾波窗口radio。

輸出:與圖像X大小相同的圖像Y。

//初始化直方圖H

//i,j分別為圖像X數組的橫縱坐標索引

//k表示濾波半徑內像素的索引

1：fori=1 toMdo

2:forj=1 toNdo

3:fork=-radiotoradiodo

4: RemoveXi+k,j-radio-1fromH;/*移除左邊一列*/

5: AddXi+k,j+radiotoH;/*增加右邊一列*/

/*新的中值在直方圖的左側，其中sum為上一次統計直方圖元素個數的和,T表示閾值*/

//判斷新的中值在直方圖的左側

6:sum>T

/*median()函數對處理好的直方圖進行尋找中值的操作*/

7:Yi,j←median(H);

/*判斷新的中值在直方圖的右側*/

8:sum

9:Yi,j←median(H);

/*判斷新的中值是否位置不變*/

10:sum=T

11:Yi,j←median(H);

12:endfor

13:endfor

14：endfor

5 實驗結果與分析

本文討論的所有算法都是在華為鯤鵬920服務器實驗環境上運行的，具體實驗環境配置如表1所示。

Table 1 Experimental environment configuration

使用C語言與ARM NEON指令集混合編程，具有良好的跨平臺性和可移植性，以下對HMPP庫的中值濾波與OpenCV的中值濾波進行性能比對分析。

5.1 8u類型下不同窗口的中值濾波結果與性能分析

本文對HMPP庫中不同的數據類型及不同的濾波窗口采用了自適應的中值濾波優化算法，其中針對小濾波窗口(3×3,5×5)的中值濾波算法優化適用所有的數據類型。圖7是在8u數據類型不同像素規模下對各算法性能進行比較的結果。后面所有性能圖中，Hmpp代表本文使用的算法，OpenCV代表對比庫OpenCV中相對應的算法。

圖7的結果表明，對于本文提出的自適應算法，對于3×3的小濾波半徑窗口，算法的最高性能達到OpenCV中值濾波算法的257%，最低為104%，對于1024×1024這種大規模圖像，性能提升達到了140%；同時對于該算法來說，除了可以處理灰度圖像，還可以處理其他數據類型的圖像，應用范圍廣，具有較高的實用價值。

Figure 7 Performance diagram of different filtering window sizes

該算法最大的創新點在于使用ARM NEON指令集進行優化，且使用了SIMD指令，算法性能提升幅度較大，達到了預期效果。對于5×5濾波窗口，HMPP庫使用的中值濾波算法的最高性能相比于OpenCV的提升了377%，最低為115%，相比于OpenCV有較大的性能提升，對于1024×1024這種大規模圖像性能提升達到了170%。

實驗結果表明，在5×5濾波窗口下，使用該算法及ARM NEON指令集優化后，相對于小濾波窗口的中值濾波，對于OpenCV具有明顯的性能優勢，整體性能大約為OpenCV的140%。

對于7×7濾波窗口，本文提出的算法性能最高達到了OpenCV的115%，最低達到了114%，針對1024×1024這種大規模圖像，性能達到了115%。可以看出，改進的直方圖性能提升明顯，同時性能不會隨著圖像規模的增大有明顯的波動，性能提升基本維持在115%。實驗結果表明，使用該算法能帶來穩定的性能提升。

對于11×11濾波窗口，從圖7中可以看到，本文提出的算法性能最高達到OpenCV的130%，最低達到了110%，針對1024×1024這種大規模圖像，性能達到了110%。實驗結果表明，在濾波窗口不斷增大時，該算法對性能不會有太大的影響。可以看出，基于直方圖的中值濾波是針對灰度圖像大濾波窗口半徑(7×7以上)或者不規則濾波窗口使用的，該算法的性能比較穩定，同時比較高效。

圖8的結果表明，當濾波窗口為3×3時，在該算法處理4通道圖像時，HMPP的中值濾波算法性能最高達到了OpenCV中值濾波算法的146%，最低為106%，對于1024×1024這種大規模圖像，性能提升達到了116%。

Figure 8 4-channel performance diagram of different filtering window sizes

實驗結果表明，HMPP多通道中值濾波算法性能穩定，對于5×5濾波窗口，HMPP庫使用的中值濾波算法的最高性能相比于OpenCV提升了140%，最低為101%，相比于OpenCV有較大的性能提升。對于1024×1024這種大規模圖像性能提升達到了101%。實驗結果表明，在5×5濾波窗口下，使用該算法及ARM NEON指令集優化后，對于小濾波窗口的中值濾波多通道，相對于OpenCV具有明顯的性能優勢，整體性能大約為OpenCV的120%。

對于7×7濾波窗口的4通道，從圖8可以看到，本文提出的算法性能最高達到了OpenCV性能的132%，最低達到130%；針對1024×1024這種大規模圖像，性能達到了130%。實驗結果表明，在多通道情況下，該算法的整體性能依舊穩定。

對于11×11濾波窗口下的4通道，從圖8可以看出，本文提出的算法性能最高達到了OpenCV性能的151%，最低達到了125%，針對1024×1024這種大規模圖像，性能達到了125%。實驗結果表明，在大濾波窗口多通道使用該算法的性能優異，相比于單通道的性能提升，多通道所帶來的性能提升幅度更大，達到了預期要求。

5.2 面向非8u中值濾波算法性能分析

本節以16u數據類型為例，分別測試不同算法在相同規模下的性能。OpenCV并沒有實現8u數據類型以上的大濾波窗口下的中值濾波算法，而本文提出的算法可以有效地解決這個問題，并提升了性能。該算法通過設置哨兵，使得隊中最小、隊尾最大，縮小每次需要比較的規模的大小，遞歸調用這個過程，求出整個濾波窗口內的中值，解決了8u以上數據類型在大濾波窗口的中值濾波問題，如圖9所示，濾波窗口大小為3×3時，本文提出的算法最高性能為OpenCV的201%，最低為150%；同樣在濾波窗口為5×5時，本文提出的算法最高性能為OpenCV的140%，最低為120%，實驗結果表明，本文采取的自適應算法整體性能優于OpenCV的中值濾波算法性能。

Figure 9 Performance of the proposed algorithm when the data type is 16u and the filtering window is 3×3 or 5×5

如圖10所示，針對大規模的圖像，本文算法的時間開銷是穩定的，也就是說，隨著規模的擴大，本文算法的優勢會明顯地體現出來，對性能有較大的提升。

Figure 10 Performance of the proposed algorithm when the data type is 16u and the filtering window is 11×11

5.3 不同中值濾波算法在不同濾波窗口大小時的性能

對于本文采取的自適應算法對中值濾波性能的提升，以8u和16u數據類型為例，在不同濾波窗口大小及不同的通道數下，進行驗證分析。本文使用8u和16類型下的3×3，5×5，7×7，11×11的性能數據。

圖11為相同大小濾波窗口下的性能差異，當濾波窗口為3×3,5×5時，本文的自適應算法性能最佳。在7×7以上的濾波窗口下，相比于歸并排序，本文使用直方圖后性能最佳，從圖11中可以看出，在不同大小的濾波窗口下，本文使用的自適應算法都能獲得較佳的性能。如圖12所示，在16u的不同大小的濾波窗口下，在小窗口下的算法相比于快速中值算法有明顯優勢，在大規模窗口下，本文采取的算法相比于歸并排序有明顯優勢。實驗結果表明，本文的自適應算法是有效的。

Figure 11 Performance comparison of different algorithms with the same filtering window when the data type is 8u

Figure 12 Performance comparison of different algorithms with the same filtering window when the data type is 16u

6 結束語

本文針對中值濾波提出了一種自適應的算法，同時使用ARM NEON指令集優化，針對不同類型的圖像和濾波窗口規模，對于8u的灰度圖像，本文提出了基于改進后的直方圖中值濾波算法，該算法相比于OpenCV中值濾波算法，性能提升了大約20%。

在濾波窗口為3×3和5×5時，使用提出的行列有序算法進行單獨處理，該算法適用于所有數據類型。

對于非8u圖像，當濾波窗口不規則或者min(H,W)>5時，本文提出了一種遞歸的算法，降低了求中值算法的時間復雜度，能穩定提升性能。

由于針對非8u大規模濾波窗口的現有算法并沒有達到很好的性能，因此后續將重點研究非8u大規模濾波窗口的高性能算法，以提升其性能。