[數字圖像處理]圖像去噪初步(1)--均值濾波器

1.圖像去噪的前言

    上一篇博文中，我對噪聲的類型進行了介紹，也使用的Matlab對各種噪聲進行了實作。舊話重提，一幅圖像，甚至是一個信号的老化，可以使用以下模型來表示。

可以使用以下算式來表示

這裡，由于退化函數

的作用，使得原圖像

産生退化(比如，運動模糊)，然後在加上一個加性噪聲項

。

      本博文，主要對去除加性噪聲的線性濾波器的性能進行了比較。對于退化函數的去除(稱為去卷積或者逆濾波)，将放在稍後的博文。

      1.1 實驗用圖像 

      1.2 實驗結果的評價

            實驗的步驟為，将實驗用圖像加上加性噪聲，然後使用濾波器進行去噪，比較所得到的圖像的畫質。這裡，就涉及到畫質的評價方法。一般的，去噪圖像的評價一般使用PSNR（峰值信噪比）。

       對于8-bit的圖檔而言，這裡的MAX為255。PSNR越大，其畫質就越好。但是，有些時候，使用PSNR來進行評價，也有不太合理的時候。

       請對比如下三張圖檔，a)是使用平均濾波器進行了處理，使其有些模糊；b)是使用高斯噪聲污染原圖；c)是使用椒鹽噪聲污染的圖像。

    問題來了，這三張圖像哪張畫質最好，哪張最差。普遍的，畫質從好到差排列，大家的答案應該是

a) > c) > b)

這樣的（更多實際例子，請參考https://ece.uwaterloo.ca/~z70wang/research/ssim/）。那麼，我們求其的PSNR是這樣的。

這明顯不科學，三幅圖像的PSNR是一樣的。反觀PSNR的計算式，PSNR計算的時候，使用了MSE這個量。而MSE僅僅表現了兩幅圖像的灰階值的差，而對于圖像的結構，卻沒有進行任何分析。

       這裡使用一種比較好的圖像畫質評價的方法：SSIM（念做：艾斯-希姆）。這是一種由兩張圖像的灰階差異，構造差異和對比度去判斷兩張圖的接近程度的方法。詳情請參考[文獻1]，這裡隻做簡單的介紹一下啦。

SSIM從圖像亮度(Luminance)，圖像對比度(Contrast)和圖像構造(Structure)去判斷處理過的圖像與原圖的差異。這裡，使用了某個區域的内的平均值作為亮度度量，使用方差作為對比度度量，使用協方差作為構造度量，來進行判斷。這樣，SSIM就比僅使用灰階去判斷的PSNR更加準确。一樣的，使用SSIM求取上面三幅圖象的類似度。

從上表可以看出來，通過使用SSIM進行判斷的結果，更加符合人眼的主觀感受。本文餘下的實驗，全部使用SSIM去判斷畫質。

2.幾個均值濾波器---線性處理

    2.1 算術均值濾波器

    算術均值濾波器很簡單，就是求某個區域内的所有像素的算術均值，可用下式子表示。     

從式子上可以看出來，這就是一個低通濾波器，會使得畫面模糊，有些許去噪能力。稍微做個實驗看看。

    将實驗用圖像加噪，噪聲均值為0、方差為0.0298的噪聲。

    下面，使用算術均值濾波器，看看去噪效果。

   被去掉了些許，隻是些許。再看頻率域内的圖像，果然是一個低通濾波器，我們都可以腦補出這個濾波器的振幅特性了，對吧？

    2.2 幾何均值濾波器

    接下來是幾何均值濾波器，求某個區域内的幾何平均值。     

對于這個濾波器，書(《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods)上說了，這個濾波器産生的去噪效果可以與算術平均值濾波器基本一樣，但是可以更少的丢失細節。     同樣的，将實驗用圖像加噪，噪聲均值為0、方差為0.17的噪聲。

       有結果可見，其去噪效果也不是太理想，但是原本晶片的pin腳什麼的，比算術均值濾波器稍微好要一些的。當然，幾何濾波器有一個緻命的缺點，一旦有0值出現，那麼這個像素的值立即被決定為0，這也就意味着，幾何濾波器不可以去除胡椒噪聲。

    2.3 算術均值濾波器與幾何均值濾波器的比較

        為了對比幾何均值濾波器與算術均值濾波器，我們進行了如下幾組實驗。由于篇幅問題，我就不貼出圖來了。看資料就可以了。         1.  噪聲：高斯噪聲，均值0，方差為0.17

         2.噪聲：高斯噪聲，均值0.2，方差為0.17

         3.噪聲：椒鹽噪聲，胡椒密度0，鹽粒密度0.1

這個還是想把結果貼上，幾何平均濾波器的實驗結果還是具有一定的觀賞性的。

         4.噪聲：椒鹽噪聲，胡椒密度0.1，鹽粒密度0

       就如同實驗結果一樣，由于包含了大量了胡椒噪聲，幾何濾波器壞掉了。得到的結果很糟糕。

       實驗結論：實驗4的資料說明 ，由于包含了0值，幾何濾波器去噪效果并不太好。但是實驗3，幾何濾波器的去噪效果真的是很不錯。簡單而言，算術均值濾波器泛用性比較好，而幾何濾波器則擅長于去除鹽粒噪聲。 ================        拓展運用：既然幾何濾波器對于去除鹽粒噪聲，那麼，對于僅僅函數胡椒噪聲的圖像取反，将胡椒噪聲轉換為鹽粒噪聲去處理，所得結果再傳回來，那麼，幾何均值濾波器還是可以去除胡椒噪聲的。 ================

    2.4 諧波均值濾波器

    其表達式如下所示，

注意式子的分母，這個濾波器不但不能去除椒鹽噪聲，對于灰階過黑的圖像而言，也是要壞事的。書(《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods)上又說了，這個濾波器，擅長于高斯噪聲的去噪。實驗一下看看。     把實驗用圖加噪，類型為高斯噪聲，均值為0，方差為0.15.

從實驗結果看，其實這個濾波器的去噪效果還不如算術平均濾波器和幾何平均濾波器。

    2.5 逆諧波均值濾波器

    逆諧波濾波器可以對應多種噪聲，式子如下。

這個濾波器，可以通過Q值的變化，來獲得一定的效果。當Q為正，這個濾波器可以去除胡椒噪聲，當Q為負，這個濾波器可以去除鹽粒噪聲。貼個結果收活吧。     

3.總結

   本文介紹了幾種圖像去噪的線性濾波器，并對他們進行了比較。本次實驗所用到的代碼，我貼在下面。    另外，評價函數ssim_index()，請自行網上查找。

close all;
clear all;
clc;

f = imread('./ckt-board-orig.tif');
f = mat2gray(f,[0 255]);
[M,N] = size(f);
%% ---------------gaussian noise-------------------
a = 0;
b = 0.1;
n_gaussian = a + b .* randn(M,N);

g_gaussian = f + n_gaussian; 
g_gaussian(find(g_gaussian > 1)) = 1;
g_gaussian(find(g_gaussian < 0)) = 0;

figure();
subplot(1,2,1);
imshow(g_gaussian,[0 1]);
xlabel('a).Image corrupted by Gaussian noise');

subplot(1,2,2);
x = linspace(-0.2,1.2,358);
h = hist(g_gaussian,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.014]),grid;
xlabel('b).The Histogram of a');
ylabel('Number of pixels');

g_diff = abs(g_gaussian - f);

MSE = sum(sum(g_diff .^2))/(M*N);
PSNR = 10*log10((1*1)/MSE)

[SSIM_G mssim] = ssim_index(f,g_gaussian,[0.01 0.03],ones(8),1);

%% ---------------Denoise(Gaussian noise)-------------------
g_Ex = zeros(M+2,N+2);
for x = 1:M
    g_Ex(x+1,:) = [g_gaussian(x,1) g_gaussian(x,:) g_gaussian(x,N)];
end
g_Ex(1,:) = g_Ex(2,:);
g_Ex(M+2,:) = g_Ex(M+1,:);

% Arithemtic mean filter
g_amf = zeros(M,N);
for x = 2:M+1
    for y = 2:N+1
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x  ,y);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x-1,y);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x  ,y-1);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x+1,y);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(  x,y+1);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x-1,y+1);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x+1,y-1);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x+1,y+1);
        g_amf(x-1,y-1) = g_amf(x-1,y-1) +  g_Ex(x-1,y-1);
        
        g_amf(x-1,y-1) = g_amf(x-1,y-1)/9;
    end
end

g_amf_diff = abs(g_amf - f);

MSE_amf = sum(sum(g_amf_diff .^2))/(M*N);
PSNR_amf = 10*log10((1*1)/MSE_amf)

[SSIM_amf mssim] = ssim_index(f,g_amf ,[0.01 0.03],ones(8),1);

% Geometric mean filter

g_gmf = zeros(M,N);
for x = 2:M+1
    for y = 2:N+1
        g_gmf(x-1,y-1) = g_Ex(x  ,y);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x-1,y);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x  ,y-1);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x+1,y);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(  x,y+1);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x-1,y+1);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x+1,y-1);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x+1,y+1);
        g_gmf(x-1,y-1) = g_gmf(x-1,y-1) *  g_Ex(x-1,y-1);
    end
end

g_gmf = (g_gmf).^(1/9);

g_gmf_diff = abs(g_gmf - f);
MSE_gmf = sum(sum(g_gmf_diff .^2))/(M*N);
PSNR_gmf = 10*log10((1*1)/MSE_gmf)

[SSIM_gmf mssim] = ssim_index(f,g_gmf ,[0.01 0.03],ones(8),1);

figure();
subplot(1,2,1);
imshow(g_amf,[0 1]);
xlabel('a).Ruselt of Denoise by Amf');

figure();
subplot(1,2,2);
imshow(g_gmf,[0 1]);
xlabel('a).Ruselt of Denoise by Gmf');

% Harmonic mean filter
g_hmf = zeros(M,N);
for x = 2:M+1
    for y = 2:N+1
        g_hmf(x-1,y-1) = 1/g_Ex(x  ,y);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x-1,y);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x  ,y-1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x+1,y);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(  x,y+1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x-1,y+1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x+1,y-1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x+1,y+1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x-1,y-1);
        
        g_hmf(x-1,y-1) = 9/g_hmf(x-1,y-1);

    end
end

g_hmf_diff = abs(g_hmf - f);
MSE_hmf = sum(sum(g_hmf_diff .^2))/(M*N);
PSNR_hmf = 10*log10((1*1)/MSE_hmf)

[SSIM_hmf mssim] = ssim_index(f,g_hmf ,[0.01 0.03],ones(8),1);

figure();
subplot(1,2,1);
imshow(g_hmf,[0 1]);
xlabel('c).Ruselt of Denoise by Hmf');

%% ---------------Denoise(salt and pepper noise)----------
close all;
clear all;
clc;

f = imread('./ckt-board-orig.tif');
f = mat2gray(f,[0 255]);
[M,N] = size(f);
%% ---------------Salt & pepper-------------------
b = 0;  %salt
a = 0.2;    %pepper
x = rand(M,N);

g_sp = zeros(M,N);
g_sp = f;

g_sp(find(x<=a)) = 0;
g_sp(find(x > a & x<(a+b))) = 1;

g_diff = abs(g_sp - f);
MSE = sum(sum(g_diff .^2))/(M*N);
PSNR = 10*log10((1*1)/MSE)

figure();
subplot(1,2,1);
imshow(g_sp,[0 1]);
xlabel('a).Image corrupted by salt&pepper noise');

%% ---------------Denoise(Salt & pepper)-------------------
g_Ex = zeros(M+2,N+2);
for x = 1:M
    g_Ex(x+1,:) = [g_sp(x,1) g_sp(x,:) g_sp(x,N)];
end
g_Ex(1,:) = g_Ex(2,:);
g_Ex(M+2,:) = g_Ex(M+1,:);

% Harmonic mean filter
g_hmf = zeros(M,N);
for x = 2:M+1
    for y = 2:N+1
        g_hmf(x-1,y-1) = 1/g_Ex(x  ,y);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x-1,y);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x  ,y-1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x+1,y);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(  x,y+1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x-1,y+1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x+1,y-1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x+1,y+1);
        g_hmf(x-1,y-1) = g_hmf(x-1,y-1) +  1/g_Ex(x-1,y-1);
        
        g_hmf(x-1,y-1) = 9/g_hmf(x-1,y-1);
    end
end

g_hmf_diff = abs(g_hmf - f);
MSE_hmf = sum(sum(g_hmf_diff .^2))/(M*N);
PSNR_hmf = 10*log10((1*1)/MSE_hmf)

figure();
subplot(1,2,1);
imshow(g_sp,[0 1]);
xlabel('a).Image corrupted by pepper noise');

subplot(1,2,2);
imshow(g_hmf,[0 1]);
xlabel('b).Ruselt of Denoise by Hmf');

% Contraharmonic mean filter
Q = -1.5;

g_cmf = zeros(M,N);
for x = 2:M+1
    for y = 2:N+1
        g_cmf_M = (g_Ex(x  ,y))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x-1,y))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x  ,y-1))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x+1,y))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(  x,y+1))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x-1,y+1))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x+1,y-1))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x+1,y+1))^(1+Q);
        g_cmf_M = g_cmf_M +  (g_Ex(x-1,y-1))^(1+Q);
        
        g_cmf_D = (g_Ex(x  ,y))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x-1,y))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x  ,y-1))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x+1,y))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(  x,y+1))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x-1,y+1))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x+1,y-1))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x+1,y+1))^(Q);
        g_cmf_D = g_cmf_D +  (g_Ex(x-1,y-1))^(Q);
        
        g_cmf(x-1,y-1) = g_cmf_M/g_cmf_D;
    end
end

g_cmf_diff = abs(g_cmf - f);
MSE_cmf = sum(sum(g_cmf_diff .^2))/(M*N);
PSNR_cmf = 10*log10((1*1)/MSE_cmf)

figure();
subplot(1,2,1);
imshow(g_cmf,[0 1]);
xlabel('b).Ruselt of Denoise by Cmf(Q=-1.5)');
           

參考文獻

[1] Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, "Image quality assessment: From error visibility to structural similarity," IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600-612, Apr. 2004.

原文發于部落格：http://blog.csdn.net/thnh169/ 

=============更新日志===================

2014.7.17     修正了一些文法錯誤。