高效均值濾波算法

摘 要:為了提高現有均值濾波算法的速度，減少冗余操作，提出了一種高效的均值濾波算法。首先建立一個輔助數組，結合濾波窗口在按行滑動時相鄰窗口之間的遞歸關系和按列移動時相鄰兩行對應的輔助數組之間的遞歸關系，設計出更新輔助數組元素和新窗口中心點對應均值的遞歸式，根據這兩個遞歸式，設計出了高效的均值濾波算法。理論分析表明該算法極大地減少了冗余操作，且算法效率不受窗口大小變化的影響。實驗結果證明，在濾波效果一致的前提下，本算法比現有快速均值濾波算法更高效。

關鍵詞:圖像平滑; 均值濾波; 滑動窗口; 算法設計

中圖分類號:TP391.41

文獻標志碼:A

文章編號:1001-3695(2010)02-0434-05

doi:10.3969/j.issn.1001-3695.2010.02.008

Highly efficient mean filtering algorithm

WANG Ke-jun， XIONG Xin-yan， REN Zhen

(College of Automation， Harbin Engineering University， Harbin 150001， China)

Abstract:For the purpose of improving the mean filter algorithms speed and reducing the redundancy operations， this paper proposed a high efficiency mean filter algorithm. First set up an assistant array. Then united the recursive relationship of adjacent filter windows， when the filter window slide along the row， and the recursive relationship of assistant arrays of adjacent two rows， when the filter window slided to the next row， to design two recursive formulas. One was for updating the items of assistant array and the other was for updating relative mean of new filter window’s center point. Using those two recursive formulas designed out the high efficiency mean filter algorithm. By theoretic analysis it could be known that the new algorithm had greatly reduced the redundancy operations and the efficiency of algorithm had nothing to do with the filter window’s size. The experiment results indicate that under the same filter effect， this algorithm is much faster than existing mean filter algorithm on processing speed.

Key words:image smoothing; mean filter; sliding window; arithmetic design

0 引言

均值濾波是一種基本的基于圖像局部統計信息對圖像進行濾波的方法，其應用廣泛。例如，利用均值濾波對圖像進行平滑處理，可以消除圖像中的點狀噪聲;利用它對圖像進行模糊操作，可以達到消除孔洞和連接曲線目標斷裂處的目的[1];經過適當的擴展，均值濾波還可以應用于紋理分割[2~4];利用均值濾波也可以對其他算法進行優化，以增加算法速度[5]。均值濾波算法與其他基于圖像局部統計信息的處理方法一樣，未經任何優化的均值濾波算法(AMF)，濾波窗口的尺寸對其效率有很大的影響。用M表示濾波窗口寬度，用N表示濾波窗口高度，對每個像素點來說，AMF算法的計算復雜度為O(M×N)。可以看出，隨著濾波窗口尺寸的增大，算法效率隨之降低。

為了減少濾波窗口尺寸變化對均值濾波算法效率的影響，以提高AMF算法的效率和穩定性，目前提出了多種改進算法。最為常用的是基于數字信號處理中常用的核分解[6]而提出的改進算法。因為均值濾波過程可以認為是利用均值濾波核對圖像信號進行線性卷積的過程，并且均值核可以分解成長度為M(濾波窗口寬度)的行向量和長度為N(濾波窗口高度)的列向量，利用它們依次對圖像進行行卷積和列卷積來完成均值濾波，這樣會使每個像素點的計算復雜度降為O(M+N)。 Gonzalez 等人[7]和Rakshit 等人[8]提出了更高效的算法。Gonzalez等人充分利用濾波窗口沿行滑動時相鄰窗口之間的遞歸關系，提出一種對每個像素點來說計算復雜度為O(2N+2)的改進算法(FMF);Rakshit 等人利用store-and-fetch[9]技術進一步減少了FMF算法的加法操作，提出一種對每個像素點來說計算復雜度為O(N+2)的算法(FMFT)。Rakshit等人還針對像素深度為整數級的圖像，提出了一種減少除法操作的方法。

本文對Rakshit 等人的FMFT方法作了進一步的改進，通過建立一個長度與待處理圖像寬度相等的輔助數組，存儲每列像素分組后的像素值之和;結合濾波窗口按行滑動時相鄰窗口之間的關系，設計出計算新窗口中心點對應均值的遞歸公式;再利用窗口中心點換行時，相鄰兩行對應的輔助數組之間的關系，設計出更新輔助數組元素的遞歸式;根據這兩個遞歸式，設計出了更高效的均值濾波算法(EFMF)。EFMF算法在降低每個像素計算復雜度的同時使得均值濾波算法效率與窗口大小無關。

本文對AMF、FMF、FMFT、EFMF四種均值濾波算法進行了實驗比較分析。應用這四種方法對不同大小的圖像進行處理，可以看出當圖像尺寸逐步增大時，與其他算法相比，EFMF算法時間消耗的增長率最低;分析了核窗口尺寸的變化對上述四種算法效率的影響，可以看到EFMF方法比其他方法效率要高，并且算法所需時間不隨著窗口尺寸的增加而變化。

1 初始均值濾波算法及改進算法

本文采用圖像處理中通常使用的坐標系統:圖像的左上角為坐標原點，橫坐標從左到右遞增，縱坐標從上到下遞增。為了方便，本文以二維灰度圖像為例對算法進行介紹，不考慮圖像邊界。

1.1 AMF算法

均值濾波算法是一種傳統的圖像處理算法，其思想是利用窗口中像素的平均值來代替窗口中心點的像素值[4]，其定義如下:

D(i，j)=∑(m，n)∈Ri，jS(m，n)/NK(1)

其中:i=0，…，H-1， j=0，…，W-1;S(i，j)表示圖像中像素點p(i，j)的灰度值;W表示圖像的寬度;H為圖像的高度;Ri，j表示均值濾波窗口(以下簡稱窗口)的中心點與圖像點p(i，j)重合時，窗口內所有圖像像素點的坐標集合，如圖1所示;D(i，j)表示目標圖像中點p(i，j)的灰度值;NK = M×N表示窗口中像素的總數。根據式(1)，初始均值濾波算法實現過程如下:

循環1:i從0到H-1

循環2:j從0到W-1

sum=∑(m，n)∈Ri，jS(m，n)

D(i，j)=sum/NK

終止2

終止1

從實現過程可以分析出，每個像素點需要M×N次加法操作和一次除法操作。要處理一幅尺寸為W×H的圖像，共需W×H×M×N次加法操作和W×H次除法操作。每個像素的計算復雜度為O(M×N)。

1.2 FMF算法

為了提高AFM算法的效率，Gonzalez等人從減少冗余操作的角度出發提出了FMF算法。該算法考慮了當核窗口中心從像素p(i，j)按行滑動到像素點p(i，j+1)時，它們對應的點集Ri，j和Ri，j+1有很大的交集Ri，j∩Ri，j+1，如圖2所示。那么在求p(i，j+1)的均值時，交集中像素點灰度值之和不用重新計算，可以通過sumj =∑(m，n)∈Ri，jS(m，n)減去Ri，j中最左邊一列像素點PLn(n=0，…，N-1)灰度值之和后得到，然后用得到的差值加上Ri，j+1中最右邊一列像素點PRn(n=0，…，N-1)灰度值之和就可以得到sumj+1 =∑(m，n)∈Ri，j+1S(m，n)。關系式表示如下:

sumj+1=sumj-∑N-1n=0S(PLn)+∑N-1n=0S(PRn)(2)

從式(2)可以看出，當窗口中心點按行滑動時，只需計算中心點像素處于第一列時的窗口內所有像素灰度值之和∑(m，n)∈Ri，0S(m，n)，同行中以后的∑(m，n)∈Ri，jS(m，n)均可利用式(2)遞歸得到。該算法的實現過程如下:

循環1: i從0到H-1

sum=∑(m，n)∈Ri0S(m，n)

D(i，0)=sum/NK

循環2:j從1到W-1

sum=sum-∑N-1n=0S(PLn)+∑N-1n=0S(PRn)

D(i，j)=sum/NK

終止2

終止1

從實現過程可以得出，圖像中第一列像素點需要M×N次加法操作和1次除法操作，如果認為加法操作效率和減法操作效率等同，其他的像素點只需2N+2次加法操作和1次除法操作。要處理一幅大小為W×H的圖像，共需H×M×N+(2N+2)×(W-1)×H次加法操作和W×H次除法操作。其中加法操作的次數可以用W×H×(2N+2)來近似。可以看出，該算法所需的加法次數比初始算法有極大程度的減少，每個像素的計算復雜度可認為是O(2N+2)，但是從中可以看出窗口高度的變化依然影響著算法的效率。

1.3 FMFT算法

Rakshit等人利用store-and-fetch技術進一步減少了FMF算法的加法操作，提出了FMFT算法。具體的方法是設計了一個長度與窗口寬度相等的先進先出(FIFO)線性隊列s，s中的每個元素存儲的是Ri，j內每列像素值之和。如圖3所示，圖中s[j]表示隊列中的第j(j=1，…，M)個元素。

FMFT算法中，對隊列s主要完成兩個操作:初始化隊列、更新隊列s。用ry表示窗口y軸方向的半徑，其初始化操作過程如下:

循環1:j從0到M-1

sum= 0

循環2 :i從-ry到ry

sum = sum+S(i，j)

終止2

將sum放入隊列s

終止1

從上面的過程可以看出，完成一次隊列初始化需要M×N次加法操作。從Ri，j到Ri，j+1時隊列s的更新過程是:s的第一個元素出列(該值即為Ri，j中最左邊一列像素值之和)，計算Ri，j+1中最右邊一列像素值之和∑N-1n=0S(PRn)放入隊列中。從此過程可知，更新一次s需要N次加法操作。

設計出數組s后，Ri，j到Ri，j+1的遞歸關系如圖4所示，計算式(2)可以變為

sumj+1=sumj-sj[1]+sj+1[M](3)

其中:sj[1]表示Ri，j對應隊列s的第一個元素;sj+1[M] 表示Ri，j+1對應隊列s的最后一個元素。與式(2)比較可知，式(3)所需的加法操作數減少了2N次。經過以上討論，FMFT算法的實現過程如下:

循環1: i從0到H-1

初始化隊列s

sum=∑Mn=1s[n]

D(i.0) = sum/NK

循環2:j從1到W-1

s的第一個元素出列，其值記為V

更新隊列s(需N次加法操作)

sum=sum-V+s[M]

D(i，j)=sum/NK

終止2

終止1

從實現過程可知，對第一列像素來說，共需H×(M×N+M)次加法操作和1次除法操作，其他像素點共需H×(W-1)×(N+2)次加法操作和一次除法操作。該算法總共需要的加法次數是H×(W-1)×(N+2)+ H×(M×N+M)= W×H×(N+2)+ H×(M×N+M-N-2)。其中H×(M×N+M-N-2)相對于W×H×(N+2)來說可忽略，也就是說該算法大約需要W×H×(N+2)次加法操作和W×H次除法操作。可以看出，該算法所需的加法次數比FMF算法有進一步的減少，每個像素的計算復雜度可認為是O(N+2)，從中可以看出FMFT相對于FMF算法，將窗口高度對算法效率的影響從2倍的關系降到了1倍，但是它還是沒有解決窗口尺寸對算法效率影響的問題。

針對像素深度為整數級的圖像，Rakshit等人也提出了一種減少除法操作的方法。下面以像素深度為256級的單通道圖像為例來介紹這種減少除法操作的方法。對整數級的圖像進行均值濾波，求得的像素點的均值也應轉換為整數，對像素深度為256級的單通道圖像而言，最后求得的像素點對應的均值meanij應為區間[0，255]內的整數，可表示如下:

meanij=round(sumRij/NK)(4)

其中:round()表示四舍五入后取整數;sumRij的取值為區間[0，255×NK]內的整數，在窗口內的所有像素點灰度值均為0時，sumRij為零，當所有的像素點灰度值為255時，sumRij為255×NK。從式(4)可以得出，只有當sumRij的值增加NK/2時，meanij的值增加1，這說明有很多不同值的sumRij對應著相同的均值。而傳統的方法是對它們都要進行一次除法操作，這樣導致了很多冗余的除法操作存在。Rakshit等人通過建立一個均值搜索數組mean，其長度為255×NK，數組mean的索引值對應著sumRij的值，數組元素中存儲的是其索引值對應的均值。建立了這個數組之后，那么每個像素對應的除法就可以避免了。數組mean的建立過程如下:

令i=j=0

循環1: n從0到NK

判段i是否等于NK/2的整數部分

是:i=0; j=j+1

判斷結束

mean[n]=j

i=i+1

終止1

通過建立均值搜索數組mean后，每個像素點對應的除法sumRi，j/NK確實可以避免，尤其是需要批處理像素深度相同的圖像時，因為均值搜索數組只與像素深度有關，所以可以共享同一個均值搜索數組。但是減少除法操作的方法帶來兩個問題:第一個問題是限制了FMFT算法適用范圍，即只能處理像素深度為整數級的圖像，而有時也需要處理像素深度為浮點類型的圖像，那么這個減少除法操作的方法就不能適用了;第二個問題是其需要很大的額外存儲空間，如當處理像素深度為16 bit的單通道圖像時，其數組長度為(216–1)×NK+1，隨著像素深度位數的增加，數組長度會呈2的指數級增大，當數組長度很大時，構建均值搜索數組mean所帶來的額外操作也不能再忽略。

2 EFMF算法

從以上的介紹可知，當前窗口內像素灰度值的和只需減去其最左邊一列的像素PLn(n=0，…，N-1)灰度值的和再加上一列新的像素PRn(n=0，…，N-1)灰度值的和，就得到了下一個窗口內像素灰度值的和。從式(2)可以看出，FMF算法采取的方法是:對每個像素點都要將這兩列的像素灰度值進行累加，從而帶來了2N次加法操作。FMFT算法通過引入一個與窗口寬度等長的存儲隊列s，從式(3)可以看到，這2N次加法操作似乎已經避免，而實際上對每個像素點它都需要對隊列s進行更新。更新一次需要N次加法操作，所以FMFT算法其實只是避免了N次加法操作，究其原因，對每個像素點來說，隊列s其實只是存儲了Ri，j最左邊一列的像素PLn(n=0，…，N-1)灰度值的和，而Ri，j+1最右邊的一列PRn(n=0，…，N-1)灰度值的和需要累加計算，從而帶來了N次加法操作。如果能先建立一個數組，其中已經分別存儲了這兩列像素灰度值的和，那么這2N次加法操作就是完全可以避免的。本文就是從這個角度對FMFT算法減少加法操作的方法作了進一步的改進。

本文提出的算法以核窗口高度N對每列像素進行分組，并求取它們像素灰度值的和存儲到數組A中，數組長度為圖像的寬度W，如圖5所示。以A[j]為例進一步說明數組s元素的意義:當核窗口中心在第i行滑動時，A[j]表示第j列中，中心點為p(i， j)、長度為N的一組像素點灰度值之和;當核窗口中心在第i+1行滑動時，A[j]表示第j列中，中心點為p(i+1， j)、長度為N的一組像素灰度值之和，如圖6所示。有了這個數組之后，再根據Ri，j和Ri，j+1的關系，用rx表示核窗口x軸方向的半徑，則有像素PLn(n=0，…，N-1)灰度值的和為A[j-rx]，PRn(n=0，…，N-1)灰度值的和為A[j+1+rx]，式(2)可以重新表示為

sumj+1=sumj-A[j-rx]+A[j+1+rx](5)

通過式(2)和(5)比較可知，式(5)比(2)少了2N次加法操作。

可以看出算法的關鍵就是建立數組A，在建立A的過程中可以充分利用行與行之間的遞歸關系減少加法操作。數組A的建立分兩個階段，即初始化和更新。其初始化過程如下:

循環1:j從0到W-1

A[j]=0

循環2 :i從-ry到ry

A[j]=A[j]+S(i，j)

終止2

終止1

從初始化過程可知，初始化需要W×N次加法操作，但是在對一幅圖像進行處理時，這種初始化過程只需一次。第二階段是對數組A的更新，以A[j]為例來說明其每個元素的更新過程。假設Ai[j]對應的是當核窗口中心在第i行時A[j]的值，當核窗口中心移動到i+1時，A[j]的值用Ai+1[j]表示，根據同一列上下兩行之間的遞歸關系，如圖6所示，Ai+1[j]可以通過Ai[j]減去像素點p(i-ry，j)的灰度值再加上像素點p(i+1+ry， j)的灰度值得到。公式表示如下:

Ai+1[j]=Ai[j]-S(i-ry，j)+S(i+1+ry，j)(6)

結合上式，當核窗口中心從第i行變到i+1行時，數組A元素的更新過程如下:

循環1:j從0到W-1

A[j]=A[j]-S(i-ry，j)+S(i+1+ry，j)

終止1

從這個過程可以看出，更新一次數組A共需2W次加法操作。

根據以上的討論，本文算法的實現過程如下:

循環1:i從0到H-1

判斷i是否為0

是:初始化數組A

否:更新數組A

判斷結束

sum=∑N-1n=0A[n]

D(i，0)=sum/NK

循環2:j從1到W-1

sum=sum-A[j-1-rx]+A[j+rx]

D(i，j)=sum/NK

終止2

終止1

從實現過程可知，對第一列像素來說共需H×N次加法操作，其他像素點總共需(W-1)×N×2次加法。考慮初始化和更新數組s，本算法總共需要的加法次數是 (W-1)×H×2+H×N+W×N+W×2×(H-1)=W×H×4+H×(N-2)+W×(N-2)。其中H×(N-2)+W×(N-2)相對于W×H×4來說可忽略，也就是說本算法大約需要W×H×4次加法操作，且所需除法操作次數為W×H。可以看到，本算法所需的操作數與窗口的大小無關，每個像素點的計算復雜度為O(4)。

Rakshit等人提出的減少除法操作的方法可以很容易地應用到EFMF算法中。采用減少除法操作的方法后，除法操作數為零，顯然與窗口大小無關，所以EFMF算法與減少除法操作的方法相結合后不會改變該算法效率獨立于窗口大小的特性。

3 實驗結果及分析

本文所介紹的算法均采用C語言實現，編譯器為GCC4.2.0，平臺CPU為Celeron 2.4 GHz，內存為512 MB，處理的對象為單通道256級Lena圖像和近紅外手背靜脈圖像，如圖7所示。

3.1 均值濾波效果比較

圖8 為AMF、FMF、FMFT、EFMF四種算法采用7×7的濾波窗口對圖8(a)進行均值濾波的效果圖。從圖8(b)~(e)圖像可以看出，這四種算法具有相同的濾波效果。用(c)~(e)分別減去(b)，得到的三幅差值圖像的像素值均為零，這定量地說明這四種算法的濾波效果相同。這是因為后三種算法只是對AMF算法進行了不同程度加法和乘法操作的優化，它們與AMF算法在原理上是完全一致的，因此并不改變最后的濾波結果。

3.2 算法效率比較

在EFMF算法不采用減少除法操作的情況下，用3×3的濾波窗口對尺寸為128×128 的單通道256級Lena圖像進行處理所需時間為0.456 2030 ms，其十分之一為0.045 620 30 ms。為方便比較和描述，以下所說的時間均為相對時間，其值為與0.045 620 30 ms的比值。

3.2.1 AMF、FMF、EFMF算法效率比較

本節所討論的EFMF算法沒有進行除法操作優化。表1記錄的是AMF、FMF、EFMF用不同的濾波窗口處理不同大小的Lena圖像所需的時間。

從表1可以看出，隨著圖像尺寸的變大，EFMF的優勢更為明顯。例如對512×512的圖像采用3×3的核進行濾波時，EFMF比FMF算法快130，比AMF算法快333;對600×800的圖像進行處理時，EFMF比FMF算法快205，比AMF算法快584。隨著濾波窗口變大，EFMF算法優勢也更為顯著。例如采用3×3的核對480×640的圖像進行處理，EFMF算法比FMF快145，比AMF快364;當采用11×11的核時，EFMF算法比FMF快684，比AMF快3 915。但是從表1中也可以看到，圖像尺寸變大時，EFMF算法的效率會有微小的變化，如768×1 024和1 024×1 024的圖像，這是因為忽略的H×(N-2)+W×(N-2)次加法操作，當圖像尺寸變大以及窗口變大時，會帶來一定的影響，但是，相對于W×H×4來說還是很小，所以說，這種影響不會太大地影響算法的效率。

圖9所示為AMF、FMF、EFMF三種算法采用7×7的濾波窗口處理不同大小的圖像算法效率的變化情況，縱坐標為時間，橫坐標為處理圖像大小與128×128的比值。從圖中可以看到，隨著圖像尺寸的增大，EFMF算法效率的優勢更加明顯。圖10所示為AMF、FMF、EFMF三種算法在處理尺寸為1 024×1 024的圖像時效率隨窗口大小變化的情況，縱坐標為時間，橫坐標為正方形濾波窗口的高度值。從圖中可以清楚地看到，EFMF算法效率要優于其他兩種算法，且其效率幾乎不受窗口大小的影響。

3.2.2 EFMF 、FMFT、FMF算法效率比較

本節討論的FMFT和EFMF算法均采用了Rakshit等人提出的減少除法操作的方法對除法操作進行了優化。圖11所示為EFMF、FMFT、FMF三種算法在處理尺寸為1 024×1 024的圖像時效率隨窗口大小變化的情況，縱坐標為時間，橫坐標為正方形濾波窗口的高度值。從圖中可以清楚地看到，FMFT算法效率要優于FMF算法，EFMF算法效率要優于FMFT算法和FMF，且其效率幾乎不受窗口大小的影響。隨著窗口尺寸的變大，EFMF算法的優勢表現更加明顯。

3.2.3 除法操作對EFMF算法效率的影響

本節將討論減少除法操作對算法效率的影響，從本文第1、2章對算法的分析可以看出，如果不減少除法操作，對一幅尺寸為W×H的圖像進行濾波，所需的總除法操作次數為W×H，只與圖像的尺寸有關，所以下面討論減少除法操作對算法效率的影響程度隨圖像尺寸變化的情況。用EFMF1表示EFMF與減少除法操作方法相結合后的算法，用EFMF2表示EFMF沒有與減少除法操作方法相結合的算法。圖12所示為EFMF1、EFMF2兩種算法采用7×7的濾波窗口處理不同大小的圖像算法效率的變化情況，縱坐標為時間，橫坐標為處理圖像大小與128×128的比值。從圖中可以看出減少除法操作后，算法的效率有所提高，且隨著圖像尺寸的變大，效率提高得越多。但是總的來看，其提高程度有限，例如處理尺寸為1 024×1 024的圖像，EFMF1所需時間為582，EFMF2所需時間為661;再考慮時間基值0.045 620 30 ms，可以看出提高的程度有限。對于減少除法操作帶來的問題在第1.3節中有討論，所以在實際的應用中應該根據需要作出選擇。

4 結束語

本文提出的算法EFMF具有效率高、冗余計算少、算法效率獨立于窗口大小等優越性能，通過與文中介紹的幾種算法進行比較可以清楚看到，EFMF算法極大地提高了均值濾波算法的效率，并且隨著圖像尺寸和濾波窗口尺寸的變大，EFMF算法在效率方面的優勢更加明顯，所以非常適合對大尺寸的圖像進行均值濾波處理。本文算法思想很容易應用到其他相關圖像處理算法中，如基于圖像的局部均值和方差的圖像二值化算法，因為求局部方差的問題可以轉換成求均值的問題，所以該算法思想可以很容易地與之結合，降低冗余計算從而提高其算法的速度。

參考文獻:

[1] PRATT W K. Digital image processing[M]. 3rd ed. New York:Wiley， 2001.

[2]PAL S K， GHOSH A， SHANKAR B U. Segmentation of remotely sensed images with fuzzy thresholding， and quantitative evaluation[J]. Remote Sensing， 2000， 21(11):2269-2300.

[3]JAIN A， ZONGKER D. Feature selection:evaluation， application， and small sample performance[J]. IEEE Trans on Pattern Analysis and Machine Intelligence，1997，19(2):153-158.

[4]章毓晉. 圖像工程:上冊[M]. 2版. 北京:清華大學出版社，2006.

[5]WANG Wen， LI Bo， ZHENG Jin， et al. A fast multi-scale retinex algorithm for color image enhancement[C]//Proc of International Conference on Wavelet Analysis and Pattern Recognition.2008:30-31.

[6]胡廣書. 數字信號處理——理論、算法與實現[M].2版. 北京:清華大學出版社，2003.

[7]GONZALEZ R C， WOODS R E. Digital image processing[M]. 2nd ed. [S.l.]:Pearson Education，2002.

[8]RAKSHIT S， GHOSH A， SHANKAR B U. Fast mean filtering technique[J]. Pattern Recognition， 2007，40(3):890-897.

[9]CORMEN H， LEISERSON C E， RIVEST R L， et al. Introduction to algorithms[M]. 2nd ed.[S.l.]:Prentice Hall of India，2001.