1.非線性處理與線性處理的差別
上一篇博文的内容，是關于均值濾波器的。比如說像算術均值濾波器，幾何均值濾波器。這以類型的濾波器的常常被用于剔除某些不需要的頻率成分，或者選擇需要的頻率成分，進而達到去噪的目的。這樣的濾波器，被稱為線性濾波器。
然而，還有一些特殊濾波器，他們被稱為非線性濾波器。其代表為中央值濾波器。所謂中央值濾波器，就是将一定範圍内的資料(對于圖像而言，是像素的灰階值)進行排序，選擇中央值作為濾波器的輸出。
中央值濾波器的目的并不是頻率成分的選擇，而是求root信号。關于root信号的定義，[文獻1]中是這樣給出的。
------------------------
定義：對于任意信号，使用中央值濾波器進行有限次處理，直到所得到的限号不變為止。中央值濾波器處理前後，如果這個信号列沒有任何變化，那麼将這個信号列稱為這個濾波器的root信号。
------------------------
對于中央值濾波器來講，可以完美去除被椒鹽噪聲污染的root信号。對于root信号，有如下的性質。對于以下定理的證明與詳細說明，請參考[文獻2][文獻3]。
=============
定理1：如果信号列{x(i)}是單調遞增或單調遞減的，那麼對于任何尺寸的中央值濾波器，這個信号列都是root信号。
定理2：如果信号列{x(i)}是在從i開始的m個信号之間的局部單調信号的話，那麼對于尺寸為N=2M+1(M<=m-2)的中央值濾波器，這樣的信号為這個濾波器的root信号。
=============
2.幾個典型的非線性濾波器
2.1 中央值濾波器 (Median Filter)
正如第一節所叙述的，中央值濾波器可以有效去除椒鹽噪聲，并且不會使得圖像模糊。我們将一副圖像染上椒鹽噪聲，然後使用不同尺寸的中央值濾波器去進行處理。所得的結果的畫質使用SSIM去評價。
首先，将圖像染上鹽粒密度為0.15，胡椒密度為0.15的椒鹽噪聲，使用3x3,5x5,7x7和9x9四種不同尺寸的中值濾波器去處理圖像。所得到的結果，如下表所示。
由上圖所示，可以看出，所用濾波器尺寸不同，得到的結果是不同的，尺寸較大的濾波器得到了不太好的結果。随着循環次數的增加，圖像的SSIM逐漸變化，最後趨近于穩定。下面是循環次數為30時候，4種尺寸濾波器得到的結果。
為了驗證噪聲與實驗結果的影響，将圖像染上鹽粒密度為0.35，胡椒密度為0.35的椒鹽噪聲，在進行濾波實驗，可以得到如下曲線。
實驗結果給出，尺寸較小的的中值濾波器，其收斂到root信号較為緩慢，其次，對于高噪聲的圖像，其修複結果是不太好的。
總結上述實驗，可以得到三個有用的結論。
1.中央值濾波器在處理過程中，不斷在尋找root信号。反複數回之後，其結果收斂，基本不改變。
2.對于不同的圖像，不同的噪聲，最适合本幅圖像的中央值濾波器的尺寸是不同的。并不是尺寸越小越好。
3.尺寸小的濾波器，收斂速度慢。反之，尺寸大的濾波器，收斂較快。
2.2 自适應中央值濾波器
我們要很好的使用中央值濾波器，隻有一個要點！就是找一個最适合的濾波器尺寸，使得原圖就是這個濾波器的root信号，那麼，就可以得到一個完美的結果（當然，最适合的濾波器是不存在的）。那麼，這裡就可以有個想法了，在執行濾波器操作的時候，盡可能多的嘗試各種尺寸，選擇一個最适合的，讓原圖很接近這個尺寸的濾波器的root信号。這樣的話，就可以得到一個不錯的結果了。
當然，之對于椒鹽噪聲的話，還可以加入這樣一個功能，由于噪聲的灰階值是已知的(鹽粒噪聲為1，胡椒噪聲0)，那麼隻有當中央值不是0或者1的時候，濾波器生效，否則，擴大濾波器尺寸。
歸納上述思想，可以得到如下的僞代碼。
START: S_xy = read(f(x,y),[size_S_xy]);
z_min = min(S_xy);
z_max = max(S_xy);
z_med = median(S_xy);
LOOP A : A_1 = z_med - z_min;
A_2 = z_med - z_max;
if ((A_1 >0) and (A_2 < 0))
goto LOOP B;
else size_S_xy++;
goto START;
if (size_S_xy <= size_max)
g(x,y) = z_med;
goto NEXT;
LOOP B : B_1 = f(x,y) - z_min;
B_2 = f(x,y) - z_max;
if ((B_1 >0) and (B_2 < 0))
g(x,y) = f(x,y) ;
goto NEXT;
else g(x,y) = z_med;
goto NEXT;
NEXT : x,y移動到下一個像素點
if (x > x_max)and(y>y_max)
return;
else goto START;
來看看這個濾波器的效果，将圖像加一個很強的噪聲椒鹽噪聲。 鹽粒密度為0.3，胡椒密度為0.3。同時将同樣的圖像，使用5x5的中央值濾波器。兩個濾波器的不同，還是很容易看出來的。
2.3 最大值濾波器與最小值濾波器
最大值濾波器與最小值濾波器這兩個濾波器，跟中央值濾波器有些類似。同樣是一個目的，用于尋找root信号，隻不過，最大值濾波器可以用于去除鹽粒噪聲，相反，最小值濾波器可以用于取出胡椒噪聲。同樣的，最大值與最小值濾波器不适用于去除，同時含有鹽粒噪聲與胡椒噪聲的圖檔。
由于這個很簡單，結果就不貼了。
2.4 中點濾波器(midpoint filter)
中點濾波器，糅合了均值濾波器和非線性的排序處理。這樣的濾波器對于去除高斯噪聲與均勻噪聲有很好的效果。
我們使用這個濾波器，與平均濾波器做個比較吧。其結果顯示在下圖。
圖檔上，縱軸表示畫質SSIM，橫軸表示噪聲的方差。綠線表示中點濾波器，紅線表示算術濾波器。可以看出，均勻噪聲上，中點濾波器的性能還是優于算術平均濾波器的。(當然，這可能跟選用的評價圖像有關，但是基本是這個趨勢)。
2.4 修正的阿爾法均值濾波器
從很小的時候開始，對唱歌節目就有這樣的一個映像，評委打分的時候，總是要去掉一個最高分，去掉一個最低分，然後求均值。這個是個什麼樣的原理呢？ 其實這是一個很典型的濾波器，可以更好的反應一個選手的水準。 在圖像處理領域，也有一個類似的處理方法，其表達式如下。
簡單的了解一下這個濾波器，首先，我們将濾波範圍
的
個資料進行排序，從大到下的順序去除
個資料，從小到大的順序去除
個資料。将剩下的資料計算均值。
這樣的濾波器，很擅長去除椒鹽噪聲與其他類型噪聲一起污染過的圖檔。依然是實驗，首先将圖檔染上鹽粒密度為0.1，胡椒密度為0.1的椒鹽噪聲。然後，再将這幅圖像染上高斯噪聲。在将圖檔進行去噪。得到的結果如下所示。
可以出，對于複數的噪聲，這樣的濾波器的效果是最好的。
3.自适應中值濾波器的Matlab代碼
由于本次代碼都比較簡單，主要是概念上的了解，是以，代碼我也就貼個自适應濾波器的代碼就完事吧！
%% ------------Adaptive median filter--------------------
close all;
clear all;
clc;
f = imread('./ckt-board-orig.tif');
f = mat2gray(f,[0 255]);
[M,N] = size(f);
% ----------------noise----------
b = 0.3; %salt
a = 0.3; %pepper
x = rand(M,N);
g_sp = f;
g_sp(find(x<=a)) = 0;
g_sp(find(x > a & x<(a+b))) = 1;
figure();
subplot(1,2,1);
imshow(g_sp,[0 1]);
xlabel('a).Image corrupted by S&P noise');
subplot(1,2,2);
x = linspace(-0.2,1.2,358);
h = hist(g_sp,x)/(M*N);
Histogram = zeros(358,1);
for y = 1:256
Histogram = Histogram + h(:,y);
end
bar(-0.2:1/255:1.2,Histogram);
axis([-0.2 1.2 0 0.04]),grid;
xlabel('b).The Histogram of a');
ylabel('Number of pixels');
[SSIM_g mssim] = ssim_index(f,g_sp,[0.01 0.03],ones(8),1);
%% Adaptive median filter
S_xy = zeros(3,3);
z_min = 1;
z_max = 0;
z_med = 0;
S_max = 7;
w = zeros(S_max,S_max);
w(round(S_max/2),round(S_max/2)) = 1;
g_Ex = imfilter(g_sp,w,'symmetric','full');
g_Admf = zeros(M,N);
S = 3;
for x = 4:M+3
for y = 4:N+3
% A
for S = 3:2:S_max
S_xy = zeros(S,S);
for s = 1:1:S
for t = 1:1:S
S_xy(s,t) = g_Ex(x+s-round(S/2),y+t-round(S/2));
end
end
z_min = min(S_xy(:));
z_max = max(S_xy(:));
z_med = median(S_xy(:));
A_1 = z_med - z_min;
A_2 = z_med - z_max;
if ((A_1>0)&&(A_2<0))
% B
B_1 = g_Ex(x,y) - z_min;
B_2 = g_Ex(x,y) - z_max;
if((B_1>0)&&(B_2<0)) g_Admf(x-3,y-3) = z_med;break;
else g_Admf(x-3,y-3) = z_med; break;
end
end
if(S >= S_max) g_Admf(x-3,y-3) = z_med; break;
end
end
end
end
figure();
subplot(1,2,1);
imshow(g_Admf,[0 1]);
xlabel('c).Ruselt of denoise by Adaptive median filter');
subplot(1,2,2);
x = linspace(-0.2,1.2,358);
h = hist(g_Admf,x)/(M*N);
Histogram = zeros(358,1);
for y = 1:256
Histogram = Histogram + h(:,y);
end
bar(-0.2:1/255:1.2,Histogram);
axis([-0.2 1.2 0 0.04]),grid;
xlabel('d).The Histogram of a');
ylabel('Number of pixels');
[SSIM_g_Admf mssim] = ssim_index(f,g_Admf,[0.01 0.03],ones(8),1);
下一篇博文，會說到去卷積處理。
參考文獻
[文獻1]非線形ディジタル信号処理 ,棟安実治 ・田口亮 著,システム制禦情報學會,朝倉書店,1999年03月20日
[文獻2]Tyan, S. G. "Median filtering: Deterministic properties." Two-Dimensional Digital Signal Prcessing II. Springer Berlin Heidelberg, 1981. 197-217.
[文獻3]Gallagher Jr N C, Wise G L. A theoretical analysis of the properties of median filters[J]. Acoustics, Speech and Signal Processing, IEEE Transactions on, 1981, 29(6): 1136-1141.
原文發于部落格：http://blog.csdn.net/thnh169/
=============更新日志===================
2014 - 7 - 17 修正了一些格式錯誤