（MATLAB/C/Python）快速中值濾波

• 一、中值濾波
• 二、快速中值濾波
• • • 介紹
    • 原理
    • 優化
• 三、代碼
• • • MATLAB
    • C
    • Python
• 四、測試
• 其他

by HPC_ZY

最近一個項目中需要用到中值濾波（不能調庫），但是核半徑相當大，用傳統的方法運行速度極慢。因此查閱文獻找到快速中值濾波的方法。寫成三種語言并分享給大家。

一、中值濾波

簡單說就是，就是對某點鄰域內所有像素進行排序，取序數在中間的值替代原始值。
這樣做對脈沖噪聲有良好的濾除作用，特別是在濾除噪聲的同時，能夠保護信號的邊緣，使之不被模糊。

實現方法
step 1：通過從圖像中的某個采樣窗口取出奇數個像素進行排序
step 2：用排序后的中值取代要處理的像素即可

代碼見后文。

注：本文主要講圖像處理，所有使用 “像素 ” 替代 “ 數據”。

二、快速中值濾波

介紹

中值濾波主要耗時就在于對窗口（鄰域）內像素進行排序，因此針對如何獲取中值進行優化。

最初——經典的中值濾波對窗口內所有像素進行排序，所以排序算法的選擇就很重要。
接著——但后來有人意識到 “我們只需要中間值”，沒必要將整個序列進行排序，因此出現了一些快速搜索數組中間值的方法。
然后——由于圖像的像素值為0~255的整數，1979年有人就提出了用累積直方圖尋找中值的方法。
最后——人們又在這基礎上不斷改進。

本文主要講基于直方圖的快速中值濾波方法

---

原理

實現方法
step 1：統計$N?N$窗口內像素
step 2：計算累積直方圖
step 3：找到累積直方圖值剛好超過$\frac{N?N}{2}$時對應的灰度值。
通過一個簡單的例子更直觀的說明。
這里假設一個只有8個灰度級的圖像，某個$5?5$的窗口內值如下，
$$?\begin{array}{ccccc}0& 4& 5& 5& 5\\ 3& 4& 7& 6& 5\\ 6& 7& 7& 7& 6\\ 3& 2& 1& 1& 0\\ 0& 1& 1& 0& 0\end{array}?$$
其統計直方圖如下表。對于5*5的窗口，中值的序數為13，所以參照累積直方圖能快速中值為4。
這里我們可以發現，其實并不需要算出整個累積直方圖，只要某一次計算超過半數，就可以停止了。

灰度01234567

直方圖	5	4	1	2	2	4	3	4
累積直方圖	5	9	10	12	14	18	21	25

分析
直方圖方法中，對于圖像灰度級為$L$、窗口尺寸為$N?N$，我們只需要進行以下操作：
1）進行$N?N$次統計（獲得直方圖）
2）最多進行$L?1$次累加+判斷（計算累積直方圖的同時判斷是否超過半數）
通常圖像的$L$都為256，我們可以推斷對于$3?3$的窗口，傳統排序方法速度更快；但當N越來越大，直方圖法的優勢就越來越大。

代碼見后文

---

優化

進一步思考可以發現，我們在平移窗口時，新窗口里有$N?(N?1)$的像素都是上一窗口里有的（如下圖黃色），如果重新統計就浪費時間。
所以有人提出不必重新計算直方圖，只用更新直方圖——移除（如下圖紅色）+加入（如下圖綠色）。
代碼見后文

---

三、代碼

本文首先用MATLAB講解演示：傳統排序中值法，快速直方圖中值法，改進快速直方圖法
然后直接分享C、Python的改進快速直方圖法。

MATLAB

經典方法
為了方便轉為C語言，沒有使用MATLAB的函數

% 經典中值濾波（不處理邊緣）（C格式） function out = medianfilterC(in,r)% 參數與初始化 [M,N] = size(in); out = zeros(M,N,'uint8'); L = 2*r+1; cidx = L*L/2+0.5; % 如果使用sort(),中值位置在這里 cache = zeros(L*L,1,'uint8'); % 遍歷圖片 for i = 1+r:M-rfor j = 1+r:N-r% 獲取鄰域for m = -r:rfor n = -r:rcache((m+r)*L+n+r+1) = in(i+m,j+n);endend% 選擇排序(找中值，所以沒必要全排序, 且沒必要保存過程的值）for p = 1:cidxminval = cache(p);idx = p;for q = p+1:L*Lif cache(q)<minvalminval = cache(q);idx = q;endendcache(idx) = cache(p);end% 中值out(i,j) = minval;end endend

---

直方圖方法

% 直方圖法（不處理邊緣）（C格式） function out = medianfilterCHist(in,r)% 參數與初始化 [M,N] = size(in); out = zeros(M,N,'uint8'); L = 2*r+1; cidx = L*L/2+0.5; % 如果使用sort(),中值位置在這里 % 遍歷圖片 for i = 1+r:M-rfor j = 1+r:N-r% 獲取鄰域直方圖hlist = zeros(256,1,'uint32');for m = -r:rfor n = -r:rtmp = in(i+m,j+n)+1;hlist(tmp) = hlist(tmp)+1;endend% 累積直方圖求中值hsum = 0;for k = 1:256hsum = hsum+hlist(k);if hsum>=cidxout(i,j) = k-1;breakendendend endend

---

改進（更新）直方圖方法

% 改進直方圖法（不處理邊緣）（C格式） function out = medianfilterCHistUpdata(in,r)% 參數與初始化 [M,N] = size(in); out = zeros(M,N,'uint8'); L = 2*r+1; cidx = L*L/2+0.5; % 如果使用sort(),中值位置在這里 % 遍歷圖片 for i = 1+r:M-r%% 計算第一列的hlist = zeros(256,1,'uint32');% 獲取鄰域直方圖for m = -r:rfor n = -r:rtmp = in(i+m,1+r+n)+1;hlist(tmp) = hlist(tmp)+1;endend% 累積直方圖求中值hsum = 0;for k = 1:256hsum = hsum+hlist(k);if hsum>=cidxout(i,1+r) = k-1;breakendend%% 計算后續的for j = 2+r:N-r% 更新直方圖for m = -r:r% 剔除tmp = in(i+m,j-r-1)+1;hlist(tmp+1) = hlist(tmp+1)-1;% 加入tmp = in(i+m,j+r)+1;hlist(tmp+1) = hlist(tmp+1)+1;end% 累積直方圖求中值hsum = 0;for k = 1:256hsum = hsum+hlist(k);if hsum>=cidxout(i,j) = k-1;breakendendend endend

---

自帶
當然了，MATLAB有他自帶的中值濾波，賊快

out = medfitler2(in, [N,N]);

---

C

改進（更新）直方圖方法

/*** 中值濾波 ***/ void medfilter2(BYTE* pImgOut, // （out）濾波后圖像（0~255）const BYTE* pImg, // （in）原始灰度圖像（0~255）int imgWidth, // （in）圖像寬（pixel）int imgHeight, // （in）圖像高（pixel）int nR // （in）窗口核半徑（pixel） ) {int nL = 2 * nR + 1;int num = nL*nL;int cidx = num / 2 + 1; // 中值所在的位置// 初始化統計直方圖int *hist = (int*)malloc(256 * sizeof(int));// 開始遍歷int tmp = 0;for (int i = nR; i < imgHeight - nR; i++){memset(hist, 0, 256 * sizeof(int));// 第一列for (int m = -nR; m <= nR; m++){for (int n = -nR; n <= nR; n++){tmp = pImg[(i + m)*imgWidth + (nR + n)];hist[tmp]++;}}int histsum = 0;for (int k = 0; k < 256; k++){histsum += hist[k];if (histsum >= cidx){pImgOut[i*imgWidth + nR] = k;break;}}// 計算后續for (int j = 1 + nR; j < imgWidth - nR; j++){for (int m = -nR; m <= nR; m++){tmp = pImg[(i + m)*imgWidth + (j - 1 - nR)];hist[tmp]--;tmp = pImg[(i + m)*imgWidth + (j + nR)];hist[tmp]++;}int histsum = 0;for (int k = 0; k < 256; k++){histsum += hist[k];if (histsum >= cidx){pImgOut[i*imgWidth + j] = k;break;}}}}}

---

Python

寫時候遇到一個小問題，由于本人剛入門不久，不太懂底層的原理。
用python跑for循環的時候運行速度特別慢，如果有大佬知道原因希望留言告訴我，謝謝了。
代碼還是放在這里。

def medfilter( img, r ):# 參數設置rows = img.shape[0]cols = img.shape[1] L = 2*r + 1num = L*Lcidx = num/2+0.5out = np.zeros([rows,cols]) # 循環求解for i in range(r, rows-r):# 計算第一列hist = [0]*256# 獲取直方圖tmp = img[i-r:i+r+1,0:L].flatten()for n in range(0,num):hist[tmp[n]] += 1hist[tmp[n]] += 1# 累積直方圖histsum = 0for k in range(0,256):histsum += hist[k]if histsum>=cidx:out[i,r] = kbreak# 后續計算for j in range(1+r,cols-r):for m in range(-r,r+1):tmp = img[i+m, j-1-r]hist[tmp] -= 1tmp = img[i+m, j+r]hist[tmp] += 1histsum = 0for k in range(0,256):histsum += hist[k]if histsum>=cidx:out[i,j] = kbreak return out;

---

四、測試

這里只展示matlab測試結果，

clear; close all; clc% 制造含噪圖像 im = imread('test.jpg'); gray = rgb2gray(im); in = imnoise(gray,'salt & pepper',0.02); % 椒鹽噪聲r = 4;% MATLAB自帶函數 tic out0 = medfilt2(in,[2*r+1,2*r+1]); toc% 經典排序方法 tic out2 = medianfilterC(in,r); toc% 直方圖中值法 tic out3 = medianfilterCHist(in,r); toc% 改進直方圖中值法 tic out4 = medianfilterCHistUpdata(in,r); toc

結果如下圖，由于我寫的沒有處理邊緣的像素，所以結果有黑邊。這里使用半徑為4，如果更大效果就非常明顯了。

其他

參考文獻
《A Fast Two-Dimensional Median Filtering Algorithm》
《Median Filter in Constant Time》
《Fast Median Filtering by Use of Fast Localization of Median Value》

由于本人剛入門Python不久，不太懂底層的原理。用python跑for循環的時候運行速度特別慢，如果有大佬知道原因希望留言告訴我，謝謝了。

總結

以上是生活随笔為你收集整理的（MATLAB/C/Python）快速中值滤波的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。