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 几何均值滤波器
2.3 算术均值滤波器与几何均值滤波器的比较
2.4 谐波均值滤波器
2.5 逆谐波均值滤波器
逆谐波滤波器可以对应多种噪声,式子如下。3.总结
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 修正了一些语法错误。