主要參考論文：Median Filter in Constant Time.pdf
參考代碼：http://files.cnblogs.com/Imageshop/CTMF.rar
中值濾波是一種經典的圖像操作，特別適用於椒鹽噪音的去除。同樣，他也是USM銳化（表示懷疑，我記得是高斯濾波)、順序處理、形態學操作（比如去孤點）等算法的基礎。更高級別的應用包括目標分割、語音和文字識別以及醫學圖像處理等。
然而，過多的處理時間嚴重的限制住了中值濾波器的使用。由於其算法的非線性和不可分離性普通的優化技術並不合適。最原始的步驟就是獲取圖像像素的一個列表，然后進行排序，接着取中值。在一般的情況下，該算法的復雜度是O(r2logr)，其中r為核的半徑。當像素的可能取值是個常數時，比如對於8位圖像，可以使用桶排序(Bucker sort），該排序使得算法的復雜度降低為O(r2)。但是除了小半徑的情況外，這樣的改進任然是不可接受的。
這里插一句，從我個人的認知上說，任何基於排序的中值濾波，都是無法對大半徑進行實時有效的處理的。
在<A Fast Two-Dimensional Median Filtering Algorithm>一文中，提出了基於直方圖統計的快速中值濾波，其時間復雜度為O(r)，這個想法我在我的Imageshop軟件里也曾經實踐過（那個時候可沒有看過這個論文，可見這個優化是個大眾都能想到的一步）。后來發現，在Paint.net里的中值用的也是這個。具體來講，這個算法通過計算在核內部的像素的直方圖間接得到中值。當從一個像素移動到另外一個與之相鄰（水平或垂直方向）的像素時，只需要更新一部分的信息，如圖1所示。為了更新直方圖，2r+1次加法以及2r+1次減法需要執行，而從直方圖中計算中值所需要的時間是一定的，如代碼段1所示。對於8位圖像，直方圖由256個元素組成，在平均上說，計算中值需要128次比較和127次加法。實際上，通過改變終止尋找的條件我們可以計算任何其它百分比效果（見代碼段1中的Percentile參數）。
圖1 黃氏算法，本圖半徑r=2
int StopAmount = (2 * Radius + 1) * (2*Radius +1) * Percentile; int Sum = 0; for (int I = 0; I <= 255; I++) { Sum += HistBlue(I); if (Sum >= StopAmount) // 滿足停止條件 { Value = I; // 得到中值 break; } }
代碼1： 從直方圖中計算中值（Percentile=0.5）或任意百分比值
為了使中值濾波的時間復雜性降低至線性以下，人們做出了很多努力。Gel使用了基於樹的算法將復雜度降低為O(log2r)，在同一篇論文中，他們簡要的說了一種復雜度為O(log r)的二維圖像中值算法。我們這里提出的算法也是用直方圖的計算來取代如龜速的排序。最近，Weiss使用層次直方圖技術獲取了O(log r)復雜度的效果，在他的方法中，雖然速度是快了，但是代碼復雜難懂。本文提出一種簡單而有效的算法，易於用CPU或其他硬件實現。
為更好的理解文章的算法，我們先來看看黃氏算法的不足。特別注意到該算法行與行之間沒有任何信息得到保留，而每個像素的處理至少有2r+1次加法和減法的直方圖計算，這就是其復雜度為O(r)的原因。憑直覺，我們猜想應該有某種方法使得對每個像素，直方圖只需累加一個固定的次數，從而獲得O(1)的復雜度。正如我們所看到的，通過保留行與行之間的信息，這變得可行。首先讓我們來介紹下直方圖的一些屬性上。對於不相鄰的區域A和B，有下式成立：
H(A ∪ B) = H(A) + H(B)
注意到兩個直方圖的累加是一個O(1)操作。他和直方圖的元素個數有關，而直方圖元素個數是由圖像位深決定的。有了這一點，我們就可以開發一個新的算法。
首先，對於每一列圖像，我們都為其維護一個直方圖（對於8位圖像，該直方圖有256個元素），在整個的處理過程中，這些直方圖數據都必須得到維護。每列直方圖累積了2r+1個垂直方向上相鄰像素的信息，初始的時候，這2r+1個像素是分別以第一行的每個像素為中心的。核的直方圖通過累積2r+1個相鄰的列直方圖數據獲取。其實，我們所做的就是將核直方圖分解成他對應的列直方圖的集合，在整個濾波的過程中，這些直方圖數據在兩個步驟內用恆定的時間保持最新。
考慮從某個像素向右移動一個像素的情況。對於當前行，核最右側的列直方圖首先需要更新，而此時該列的列直方圖中的數據還是以上一行對應位置那個像素為中心計算的。因此需要減去最上一個像素對應的直方圖然后加上其下面一像素的直方圖信息。這樣做的效果就是將列直方圖數據降低一行。這一步很明顯是個0(1)操作，只有一次加法和一次減法，而於半徑r無關。
第二步更新核直方圖，其是2r+1個列直方圖之和。這是通過減去最左側的列直方圖數據，然后再加上第一步所處理的那一列的列直方圖數據獲得的。這一步也是個O(1)操作，如圖2所示。如前所述，加法、減法以及計算直方圖的中值的耗時都是一些依賴於圖像位深的計算，而於濾波半徑無關。
圖 2 算法的兩步執行 (a)右側的列直方圖的更新通過減最上部和加最下部像素信息獲得
(b) 核直方圖的更新通過減最左側和加最右側列直方圖信息獲得
上述的實際效果就是核直方圖向右移動，而列直方圖向下移動。在計算中，每個像素只需訪問一次，並且被添加到一個直方圖數據中。這樣，最后一步就是計算中值了，如代碼段1所示，這也是一個O(1)操作。
綜上所述，所有的單像素操作（包括更新列以及核直方圖、計算中值）都是 O(1)操作。現在，我們重點來說說初始化操作，即通過累積前r行的數據來計算列直方圖以及從前r列直方圖數據計算第一個像素點的核直方圖。這個過程是個O(R)操作。此外，從一行移動到下一行也占用了另外一個O(R)操作(表示不理解)。然而，既然這個初始化操作對每行只執行一次，相對於其他計算來說，這個可以忽略不計。
針對8位灰度圖像，我們對上述算法進行一下總結。
（1）、對核最右側的列直方圖執行一次加法。
（2）、對同一列直方圖執行一次減法，去除多余的像素信息。
（3）、將更新后的列直方圖數據加到核直方圖中，這進行了256次加法。
（4）、將無效的列直方圖數據從核直方圖中減去，這需要256次減法。
（5）、為找到核直方圖的中值，平均需要128次比較和127次加法。
上述計算量看起來比較多。然而，大部分操作都可並行處理，從而明顯的降低了處理時間。更重要的，還有很多優化措施能降低計算量，如下所示。
1、並行化
現代處理器提供了SIMD指令可以用來加速我們的算法。上面描述的操作大部分都是對直方圖數據進行加和減操作。通過MMX,SSE2或Altivec指令可以並行處理多個直方圖操作。為了在一條指令中做更多的直方圖處理，直方圖的只能用16位的數據類型。因此，核心的大小最大為216(對應的核的大小是（256/2-1))，對於一般的應用足夠了。這個限制並不只是對我們的算法：這是一種優化的手段而已。
另外一個可以運行並行的地方就是從圖像中讀取數據以及將其累加到對應的直方圖中。同上述交替更新列和核直方圖不同的是，我們可以首先更新整行的列直方圖。然后利用SIMD指令，我們可以並行的更新多列直方圖數據（不同列直方圖的更新之間沒有任何關系，所以好並行）。然后再像平常一樣更新核直方圖。
這條優化手段對於有些高級語言是無法實現的。像VC6，VC.NET這類可以直接內嵌匯編的語言，雖然可以實現，也需要作者具有很好的匯編語言基礎，因此，實施的難度比較大。有興趣的讀者可以參考附件中的中的SSE代碼。
2、緩存優化
恆常時間的中值濾波算法需要在內存中為每列保持一個直方圖，對於圖像，這很容易就多達數百KB的大小，通常這大於今天的處理器的緩存。這導致訪問內存的效率降低。一種方式就是將圖像在水平方向上分成幾部分分開處理。每個塊的寬度大小需精心選擇， 使得其對應的列直方圖不大於緩存的大小而有能充分利用緩存。這種修改的不利點就是擴大的邊緣效應。在實踐中，這通常能導致處理時間的大幅降低。請注意，在不同的處理器上同時處理這些塊是該算法的一種很簡單的並行算法。
這種優化說實在我不知道如何用代碼去實現。
3、多層直方圖
在<A Coarse-to-Fine Algo-rithm for Fast Median Filtering of Image Data With a Huge Number of Levels>一文中顯示了多尺度直方圖是非常有效的優化手段。其想法是維持一個平行的較小的直方圖，直方圖記錄了圖像的高位數據。例如，對於8位圖像，使用兩層的直方圖很常用，其中高層使用4位，而低層使用全8位數據。習慣上我們分別給他們命名為粗分和細分直方圖。粗分直方圖包含了16(2^4)個元素,每個元素是對應的細分直方圖高位的累積和。
使用多層直方圖有兩個好處，第一個就是計算中值過程的加速。我們可以首先在粗分數據中需找到中值在細分數據中段的位置而不用檢查整個256個位置。平均上說這只需要16次而不是128次比較和相加。第二個好處是關於直方圖的相加和相減。當對應的粗分數據為0時，則可以不用計算對應的細分段。當半徑 r很小時，列直方圖是稀疏分布的，這個時候的分支判斷是很有必要的。
以上說的很籠統。舉個簡單的例子吧，對於3*3的一個核，假如像素的數據值分別為： 100、120、98、77、215、50、243、199、180。對應的高位序列為：6、7、6、4、13、3、15、12、11。低位序列為：4、8、2、13、7、2、3、7、4。 那么粗分直方圖的16個元素的值分別為:
Coarse[3]=1 Coarse[4]=1 Coarse[6]=2 Coarse[7]=1 Coarse[11]=1 Coarse[12]=1 Coarse[13]=1 Coarse[15]=1，其他都為0；
中位數的累加值為3*3/2=5，對粗分直方圖進行累加：Coarse[3]+Coarse[4]+Coarse[6]+Coarse[7]=5，得到中值對應的高位段是H=7，因此，中間值就應該在112-128之間，然后再從細分直方圖的對應段中按照類似的方式找到低位的值L，最后得到中值H*16+L。
4、條件更新核
粗分和細分直方圖的分離還有一個不明顯但是很有效的優化效果。注意到算法的大部分時間都耗費更新在核直方圖的時加上或減去列直方圖數據，這個時間隨着實時更新粗分直方圖而有條件的更新細分直方圖而得到降低。
記得前面說過計算中值的過程是先在粗分數據中尋找中值所在段，然后再從細分數據中找到精確值。對於核的中值，每個列直方圖最多只會有2r+1次貢獻，意味着只有2r+1個對應的細分段對計算結果有用。那些從來未被使用的段，其對應的細分數據將無需更新。
為了實現該功能，我們需要為每個開辟一個記錄其最后被更新的位置的列表。當從一個像素移向下個一個像素時，我們更新列直方圖以及核直方圖的粗分數據。然后根據粗分數據計算出中值再細分數據中所在的段。下一步，根據這個段上次被更新的位置更新的細分直方圖。如果上次更新的位置和當前列的位置相差2r+1的距離，那說明舊的位置和當前位置沒有任何交叉。因此我們從0開始更新段。正是通過這種方式，我們加速了整個處理速度。
交錯布置直方圖數據，從而使得相鄰列的直方圖數據在內存也是相鄰的是有好處的。因此，細分直方圖數據需要按下述方式布置：段索引、列索引、最后是元素索引。這樣，核的細分直方圖的更新就是對一塊連續的內存的累加了，具體的講，細分直方圖有類似如下的定義形式:int [,,] Fine= new int [16,Width,16],其中第一個16對應段索引，即像素值的高4位，最后一個16是元素值，即像素的低4位值。因為三維數組的訪問會造成冗余的計算下標的過程，因此，為了提高速度，應該使用一維數組或者直接用指針訪問內存，以一維數組為例，此時的定義應該修改為int [] Fine=new int[16*Width*16],對於第X行某個像素值Value，其對應的細分位置則為： Fine[(Value>>4) * (Width*16)+X*16+ (Value & 0xF)];
具體的可能還是看參考代碼更容易理解。
參考代碼中如下函數：
static inline void histogram_sub( const uint16_t x[16], uint16_t y[16] ) { int i; for ( i = 0; i < 16; ++i ) { y[i] -= x[i]; } }
第一個建議是把這個小循環手工展開。第二，我是是用C#編程實現結果的，C#沒有inline，於是我把這樣的代碼直接展開內嵌到我的代碼中，可是令人詫異的結果卻是調用函數的版本速度卻比內嵌的速度還要快，反匯編函數版本代碼如下：
y[0] -= x[0]; 00000000 push ebp 00000001 mov ebp,esp 00000003 push esi 00000004 mov esi,ecx 00000006 mov ecx,edx 00000008 mov eax,dword ptr [esi] 0000000a sub dword ptr [ecx],eax y[1] -= x[1]; 0000000c lea edx,[ecx+4] 0000000f mov eax,dword ptr [esi+4] 00000012 sub dword ptr [edx],eax y[2] -= x[2]; 00000014 lea edx,[ecx+8] 00000017 mov eax,dword ptr [esi+8] 0000001a sub dword ptr [edx],eax ...........................其他的減法..................................... y[14] -= x[14]; 00000074 lea edx,[ecx+38h] 00000077 mov eax,dword ptr [esi+38h] 0000007a sub dword ptr [edx],eax y[15] -= x[15]; 0000007c add ecx,3Ch 0000007f mov edx,ecx 00000081 mov eax,dword ptr [esi+3Ch] 00000084 sub dword ptr [edx],eax 00000086 pop esi } 00000087 pop ebp 00000088 ret
關於這個問題，我的分析是，寫成函數的版本，雖然多了幾句壓棧和出棧的語句外，CPU會充分利用寄存器來進行操作，而內嵌后，由於變量太多，CPU只能利用內存來處理這些來處理這些賦值的。而寄存器比內存的操作快多了。因此不曉得VC的內聯函數會不會也有這問題。
關於算法的耗時情況，原文給出了一個圖表：
我本機上安裝的是PS CS4，經過測試，其中間值算法的耗時也是隨着用戶輸入的半徑的增加而成線性增加的，但是在小半徑的時候還是比較快的。
對於小半徑，我的建議是采用黃算法的模式+多層直方圖的方式來實現，速度會比本文要更快些。從理論上分析，而只有在半徑大於7時，可采用本文效果達到O(1)復雜度。
原始圖像 半徑=5，百分比=50 半徑=40，百分比=50
半徑=5，百分比=25 半徑=5，百分比=75 半徑=40，百分比=75
以一副1024*768的24位真彩色圖像為例，進行速度結果匯報：
| 半徑 | 本文算法 | 黃算法+多層直方圖 |
| 2 | 506 | 259 |
| 4 | 512 | 323 |
| 6 | 478 | 377 |
| 8 | 469 | 452 |
| 10 | 479 | 515 |
| 20 | 467 | 1004 |
| 50 | 483 | 2333 |
| 100 | 525 | 4947 |
同樣，提供個編譯好的文件給有興趣研究該算法的朋友看看效果：
http://files.cnblogs.com/Imageshop/MedianFilter.zip
由於上述MedianFilter是用VB6編制的，而VB6不支持指針，只能用數組來操作圖像數據，在有些地方會造成效率的嚴重丟失，我用C#寫的速度（未用多線程）是該MedianFilter的速度的三倍到4倍左右，如果是處理灰度，基本上又可一達到同樣大小彩色圖像速度的2到3倍。
在實際應用中，大半徑的中值對於去除椒鹽噪音是沒有意義的，因為此時圖像已經損失了太多有用的信息了。根據我的了解，大半徑可以發揮用處的地方有：1、如果你的程序有和PS一樣的選區技術，那么選區的平滑這個功能其實就是對選區數據進行中值處理的過程，這個當然希望之星速度和半徑無關。2、一些二值圖像的去噪上可以利用，比如一定半徑的中值，可以去除一些孤立的塊，從而消除噪音。3、在對某些圖像進行藝術處理時，需要大半徑的中值。
2013.10.16 補充：
最近開始學習了下C，於是用C語言實現上述過程。而在論文作者的原始代碼中，用到SSE2的相關函數進行處理，即有如下的代碼：
__inline void HistgramAdd( unsigned short *x, unsigned short *y ) { *(__m128i*) y = _mm_add_epi16( *(__m128i*) y, *(__m128i*) x ); *(__m128i*) (y+8) = _mm_add_epi16( *(__m128i*) &y[8], *(__m128i*) &x[8] ); } __inline void HistgramSub(unsigned short *x, unsigned short *y ) { *(__m128i*) &y[0] = _mm_sub_epi16( *(__m128i*) &y[0], *(__m128i*) &x[0] ); *(__m128i*) &y[8] = _mm_sub_epi16( *(__m128i*) &y[8], *(__m128i*) &x[8] ); }
這里_mm_add_epi16是一組Instructions指令，編譯器在編譯的時候會將其編譯為SSE對應的匯編代碼，而由於這些指令能實現指令級別並行，比如上述_mm_add_epi16可以在同一個指令周期對8個16位數據同時進行處理，並且HistgramAdd這些函數在程序里會大量使用到，因此程序的速度能大幅提高。
對於上述測試的1024*768的圖片，用C編寫DLL，C#調用，同樣的機器耗時平均在160ms左右，速度較純粹的C#的代碼快了3倍多，可見SSE的強大。
[DllImport("MedianBlur.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Unicode, ExactSpelling = true)] private static extern void MedianBlurGray( byte * Src, byte * Dest, int Width, int Height ,int Stride ,int Radius,int Percent); [DllImport("MedianBlur.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Unicode, ExactSpelling = true)] private static extern void MedianBlurRGB( byte * Src, byte * Dest,int Width, int Height ,int Stride ,int Radius,int Percent);
附件為C#調用C編寫的DLL的工程代碼文件：http://files.cnblogs.com/Imageshop/MedianBlurTest.rar
由於_mm_add_epi16是針對16位的數據類型進行的處理，所以中值得半徑一般要求不大於128，否則數出現數據溢出等錯誤，工程中這么大的半徑已經足夠應付大部分場合的。
***************************作者： laviewpbt 時間： 2013.4.26 聯系QQ: 33184777 轉載請保留本行信息*************************