图像处理第三章作业3
题目
- 仿真实现多幅图像平均去高斯白噪声
- 仿真实现均值滤波去除高斯白噪声
- 仿真实现中值滤波去除脉冲噪声
- 分别用拉普拉斯算子和sobel算子实现图像的锐化增强，并对比实验结果
实验代码、结果及其分析
1.平均除高斯白噪声
代码
clear;
X = imread('image1.jpg');
X = rgb2gray(X);
[m,n] = size(X);
sigma = 0.09;% 产生噪声的方差
subplot(3,3,1)
imshow(uint8(X))
title('原图')
X0 = im2double(X); % 为方便计算转化成浮点数
for j = 1:8
X1 = zeros(m,n);
for i = 1:2^j
X1 = X1 + X0 + randn(m,n) * sqrt(sigma);
end
X2 = X1./(2^j);
subplot(3,3,j+1)
imshow(X2)
title(['k=',num2str(2^j)])
end
结果
分析
代码实现思路很简单，将带噪声的图片多次叠加之后除以张数，就可以得到对应结果。从结果中可以看出，k值（叠加次数）越大，去高斯白噪声的效果越好。k=256时噪声已经不太明显了。
2.均值滤波去除高斯白噪声
代码
clear;
X = imread('image1.jpg');
X = rgb2gray(X);
[m,n] = size(X);
figure(1)
subplot(3,3,1)
imshow(uint8(X))
title('原图')
X0 = im2double(X); % 为方便计算转化成浮点数
X0 = X0 + randn(m,n)*sqrt(0.09);
subplot(3,3,2)
imshow(X0)
title('加噪声后')
for k = 1:7
X1 = zeros(m+4*k,n+4*k);
X1(2*k+1:2*k+m,2*k+1:2*k+n) = X0; % 边界处理 补0
B1 = ones(4*k+1,4*k+1);
for i = 2*k+1:(2*k+m)
for j = 2*k+1:(2*k+n)
s = sum(sum(X1(i-2*k:i+2*k,j-2*k:j+2*k).*B1))/(4*k+1)^2; % mask与原矩阵相乘并得到总和的平均值
X1(i,j) = s; % 中心位置的像素替换，其他保持原图
end
end
subplot(3,3,k+2)
imshow(X1)
title(['mask:',num2str(4*k+1),'*',num2str(4*k+1)])
end
figure(2)
subplot(2,2,1)
imshow(uint8(X))
title('原图')
X0 = im2double(X); % 为方便计算转化成浮点数
X0 = X0 + randn(m,n)*sqrt(0.09);
subplot(2,2,2)
imshow(X0)
title('加噪声后')
X1 = zeros(m+2*1,n+2*1);
X1(2:m+1,2:n+1) = X0; % 边界处理 补0
B2 = [1,2,1;2,4,2;1,2,1];
B3 = [1,1,2,1,1;1,2,4,2,1;2,4,8,4,2;1,2,4,2,1;1,1,2,1,1];
for i = 2:(m+1)
for j = 2:(n+1)
s = sum(sum(X1(i-1:i+1,j-1:j+1).*B2))/sum(sum(B2)); % mask与原矩阵相乘并得到总和的平均值
X1(i,j) = s; % 中心位置的像素替换，其他保持原图
end
end
subplot(2,2,3)
imshow(X1)
title('3*3加权')
X2 = zeros(m+2*2,n+2*2);
X2(3:m+2,3:n+2) = X0; % 边界处理 补0
for i = 3:(m+2)
for j = 3:(n+2)
s = sum(sum(X2(i-2:i+2,j-2:j+2).*B3))/sum(sum(B3)); % mask与原矩阵相乘并得到总和的平均值
X2(i,j) = s; % 中心位置的像素替换，其他保持原图
end
end
subplot(2,2,4)
imshow(X2)
title('5*5加权')
结果
分析
这个实验分别对均值滤波和加权均值滤波效果都做了仿真，同时对边界做了补零处理，遍历图中所有像素，以之为中心点与均值滤波模板进行矩阵相乘，取总和均值后将处理之后的值放回到原像素位置上。可以看到高斯白噪声确实一定程度上被去除，并且随着模板越来越大，去除噪声效果也逐渐变好，但是图片会因此变模糊。对比加权均值滤波和均值滤波会发现其实两种处理对测试图片来说区别不大，但从概念中可以知道加权均值滤波应该可以在去除噪声前提下令原图没那么模糊，更忠于原图片数据。
3. 中值滤波去除脉冲噪声
代码
clear;
X = imread('image1.jpg');
X = rgb2gray(X);
[m,n] = size(X);
figure(1)
subplot(3,3,1)
imshow(uint8(X))
title('原图')
X0 = im2double(X); % 为方便计算转化成浮点数
X0 = X0 + round(randn(m,n)*sqrt(0.09));% 四舍五入将数据都变成0或1
subplot(3,3,2)
imshow(X0)
title('加噪声后')
for k = 1:7
X1 = zeros(m+4*k,n+4*k);
X1(2*k+1:2*k+m,2*k+1:2*k+n) = X0; % 边界处理 补0
for i = 2*k+1:(2*k+m)
for j = 2*k+1:(2*k+n)
A = X1(i-2*k:i+2*k,j-2*k:j+2*k);
A = reshape(A,1,(4*k+1)^2); % 把矩阵展开成一行向量
B = sort(A); % 排序
s = B(ceil((4*k+1)^2/2)); % 取中值（因为这里都是奇数，所以直接除2向上取整就是中值）
X1(i,j) = s; % 中心位置的像素替换，其他保持原图
end
end
subplot(3,3,k+2)
imshow(X1)
title(['mask:',num2str(4*k+1),'*',num2str(4*k+1)])
end
结果
分析
中值滤波中也对图像做了边界补零处理，但是从结果图可以看出，当中值滤波模板过大之后，会对图像效果造成毁灭性的打击，当中值滤波模板为5*5时，去噪声效果最好，图像也几乎与原图无异。但模板大小增大后，我们可以看到，图像灰度的非连续性变差，导致图片在另外一种意义上也变模糊了，到最后甚至变为全黑。但是笔者在未做边界处理（如图3）时，可以看到虽然中值滤波模板大小与图1中的一致，但是随着模板大小增大并不会最后令图片全黑，但是一开始的处理效果（5*5之间的对比）也会是做了边界处理之后的效果更好，笔者这里推测由于中值滤波的特性，当我们做了边界补零处理后，黑色像素点会对中值滤波效果造成较大的影响。而图2则是用均值滤波处理椒盐噪声的效果，可以看到与前一个处理高斯噪声相比效果没那么好，而与中值滤波效果对比也明显是中值滤波效果要好很多。
4.分别用拉普拉斯算子和sobel算子实现图像的锐化增强，并对比实验结果
代码
clear;
X = imread('image1.jpg');
X = rgb2gray(X);
[m,n] = size(X);
X0 = im2double(X); % 为方便计算转化成浮点数
sobelx = [-1,0,+1;-2,0,2;-1,0,1];
sobely = sobelx';
gs1 = [0,1,0;1,-4,1;0,1,0];
gs2 = -gs1;
gs3 = [1,1,1;1,-8,1;1,1,1];
gs4 = -gs3;
gs5 = [0,-1,0;-1,5,-1;0,-1,0];
gs6 = [-1,-1,-1;-1,9,-1;-1,-1,-1];
figure(1)
subplot(3,3,1)
imshow(uint8(X))
title('原图')
subplot(3,3,2)
X1 = myImFilter(X0,ones(3,3)./9);
title('均值滤波（模糊）后')
subplot(3,3,3)
A = myImFilter(X1,gs1);
title('中心-4的拉氏算子')
subplot(3,3,4)
myImFilter(X1,gs2);
title('中心+4的拉氏算子')
subplot(3,3,5)
imshow(X1-A)
title('锐化处理后效果1')
subplot(3,3,6)
B = myImFilter(X1,gs3);
title('中心-8的拉氏算子')
subplot(3,3,7)
imshow(X1-B)
title('锐化处理后效果2')
subplot(3,3,8)
myImFilter(X1,gs5);
title('中心为+5扩展锐化模板')
subplot(3,3,9)
myImFilter(X1,gs6);
title('中心为+9扩展锐化模板')
figure(2)
subplot(3,3,1)
imshow(uint8(X))
title('原图')
subplot(3,3,2)
imshow(X1)
title('均值滤波（模糊）后')
subplot(3,3,3)
AX = myImFilter(X1,sobelx);
title('X方向的sobel算子')
subplot(3,3,4)
AY = myImFilter(X1,sobely);
title('Y方向的sobel算子')
subplot(3,3,5)
imshow(X1-AX)
title('X方向处理后效果')
subplot(3,3,6)
imshow(X1-AY)
title('Y方向处理后效果')
subplot(3,3,7)
imshow(X1-AX-AY)
title('XY方向处理后效果')
%myImFilter.m
function [afterImg]=myImFilter(img,mask)
[m,n] = size(img);
[a,b] = size(mask);
T = zeros(m,n); % 保证处理前后大小一致，在外面一层补0
for i = 1:(m-a+1)
for j = 1:(n-b+1)
s = sum(sum(img(i:i+a-1,j:j+b-1).*mask)); % mask与原矩阵相乘并得到总和
T(i+2*floor(a/2),j+2*floor(b/2)) = s; % 提取边缘
end
end
afterImg = T;
imshow(afterImg)
end
结果
分析
这里为了便于用各种算子去处理图片，自己另外手写了一个模仿imfilter的函数。从对比中我们可以看到，中心为+4和-4的拉普拉斯算子提取的边缘效果是一致的，这也符合它的旋转不变性。而-8与-4的算子比较可以看出-8提取边缘细会比-4的更要多一些，锐化后的处理效果也相对于模糊处理后的图片有了明显效果，-8的锐化明显是比-4的更强。而扩展拉普拉斯算子的处理相当于直接对原图进行锐化，输出的就是一张锐化处理后的图片而非边缘提取。可以看出中心为+5的扩展和-4提取边缘加原图后的处理效果是一致的，同理+9和-8也是如此。而对比用sobel算子提取x，y方向的边缘效果，可以发现sobel算子提取边缘效果更好（图中人物的五官也提取了出来），相对应的其最终锐化处理效果也更为良好。笔者认为这也是不同算子的特性所决定的，sobel算子的分别提取x，y方向梯度后合并，确实会比直接以中心为基准的拉普拉斯算子效果会更好。
总结和感想
这次实验亲自手写仿真了多种降噪和锐化过程，从笔者以往图片处理的经验来看，降噪和锐化是非常常用的两个手段，大多数图片都需要经过这两种处理才能得到比较优美的观赏性，因此通过这次实验对这两种图像处理手法有了进一步的了解，与实际生活中的图像处理联系起来感觉有了不一样的体会。
