【数字图像处理】【Matlab】【实验二】图像增强（1）

author：zox

实验环境：Matlaba

图像增强（1）

一、实验目的二、实验题目三、实验内容3.1 相关知识1.噪声2.空域滤波3.中值滤波4.均值滤波`5.加墙+卷积操作(convolution)`6.锐化算子7. 实验中直接使用的函数 3.2 实验代码【sy2.m】【mymiddle.m】【myaverage.m】【mysharpen.m】【myconvolve.m】 3.3 实验结果四、实验心得

一、实验目的

熟悉并掌握MATLAB图像处理工具箱的使用；理解并掌握常用的图像增强技术。

二、实验题目

为图像分别加高斯噪声和椒盐噪声，采用中值滤波方法对受噪声干扰的图像滤波，窗口分别采用 3 ∗ 3 , 5 ∗ 5 , 7 ∗ 7 3*3,5*5,7*7 3∗3,5∗5,7∗7。对受噪声干扰的图像进行均值滤波。分别采用 R o b e r t s Roberts Roberts 算子、 S o b e l Sobel Sobel算子、 P r e w i t t Prewitt Prewitt算子和 L o g Log Log算子进行图像锐化。

三、实验内容

3.1 相关知识

1.噪声

数字图像在获取、传输的过程中都可能会受到噪声的污染，常见的噪声主要有高斯噪声和椒盐噪声。

其中，高斯噪声主要是由摄像机传感器元器件内部产生的；椒盐噪声主要是由图像切割所产生的黑白相间的亮暗点噪声，“椒”表示黑色噪声，“盐”表示白色噪声。

2.空域滤波

空域滤波是在图像空间中借助模板对图像进行邻域操作，处理图像每一个像素的取值都是根据模板对输入像素相应领域内的像素值进行计算得到的。

空域滤波基本上是让图像在频域空间内某个范围的分量受到抑制，同时保证其他分量不变，从而改变输出图像的频率分布，达到增强图像的目的。

空域滤波器是应用模板卷积方法对图像每一像素进行局部处理。模板就是一个滤波器，使用不同的滤波器可产生不同的滤波效果。

空域滤波分为线性滤波和非线性滤波。

3.中值滤波

中值滤波就是将窗口中奇数个数据按大小顺序排列,处于中心位置的那个数作为处理结果。

在一定条件下，克服线性滤波器（如邻域平滑滤波等）的模糊，对滤除脉冲干扰及图像扫描噪声最为有效。不需要统计特性，方便。但是对一些细节多，特别是点、线、尖顶细节多的图像不宜采用中值滤波。

具体过程：设有一个一维序列 f 1 , f 2 , … , f n f1, f2, …, fn f1,f2,…,fn，取窗口长度m（ m为奇数），对其进行中值滤波，就是从输入序列中相继抽出m个数，再将这m个点按其数值大小排序，取其序号为中心点的那个数作为滤波输出。

4.均值滤波

通过一点和邻域内像素点求平均来去除突变的像素点，优点是算法简单，计算速度快，代价是会造成图像一定程度上的模糊。

具体流程：当灰度图像f中以像素f(i， j)为中心的N×N屏蔽窗口(N=3，5， 7…)内平均灰度值为a时，无条件作f(i， j)=a处理， N由用户给定，且取N值越大，噪声减少越明显。但平均是以图像的模糊为代价的。

`5.加墙+卷积操作(convolution)`

（1）加墙

通常来说，图像的周围都是补零或者补原来的边沿。这里我们选择补充图像原来的边沿；

其次边沿应该补充多少呢？通过对我们通常所选择的模板3×3、5×5的进行分析，可得到当我们选择（k×k）[其中k为奇数]的模板时，对应的“墙”厚度应该为" k − 1 2 \frac{k-1}{2} 2k−1",下图为原图为2×2时，使用3×3的模板的补墙过程。

图1 补墙过程

然后对于任意一点，如点（1，1）与其对应模板的关系如下图所示：

图2 模板与中心点的关系

则对于（k×k）的模板，其中心点标号为" k ∗ k + 1 2 \frac{k*k+1}{2} 2k∗k+1"

由此便可通过求这k×k个点的像素值的中值和均值分别实现中值滤波和均值滤波。

（2）卷积操作

当我们选择（3×3）的模板时，卷积操作计算为：

Z = [ z 1 z 2 z 3 z 4 z 5 z 6 z 7 z 8 z 9 ] Z= \begin{bmatrix} z1 & z2 & z3 \\ z4 & z5 & z6 \\ z7 & z8 & z9 \\ \end{bmatrix} Z=⎣⎡z1z4z7z2z5z8z3z6z9⎦⎤ m o d e l = [ k 1 k 2 k 3 k 4 k 5 k 6 k 7 k 8 k 9 ] model= \begin{bmatrix} k1 & k2 & k3 \\ k4 & k5 & k6 \\ k7 & k8 & k9 \\ \end{bmatrix} model=⎣⎡k1k4k7k2k5k8k3k6k9⎦⎤

z 5 = z 1 ∗ k 1 + z 2 ∗ k 2 + z 3 ∗ k 3 + … … + z 8 ∗ k 8 + z 9 ∗ k 9 z5=z1*k1+z2*k2+z3*k3+\ldots\ldots+z8*k8+z9*k9 z5=z1∗k1+z2∗k2+z3∗k3+……+z8∗k8+z9∗k9

选择其他大小模板时操作类似，可参照加墙的操作来理解。根据模板的不同，图像卷积操作后的效果也不同。

6.锐化算子

1） R o b e r t s Roberts Roberts算子：利用局部差分算子寻找边缘，边缘定位精度较高，但容易丢失一部分边缘，同时由于图像没经过平滑处理，因此不具备能抑制噪声能力。对陡峭边缘且含噪声少的图像效果较好。

首先使用以下模板分别对图像进行卷积操作：

m o d e l x = [ 0 0 0 0 − 1 0 0 0 1 ] m o d e l y = [ 0 0 0 0 0 − 1 0 1 0 ] \ \ modelx=\left[\begin{matrix}0&0&0\\0&-1&0\\0&0&1\\\end{matrix}\right]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ modely=\left[\begin{matrix}0&0&0\\0&0&-1\\0&1&0\\\end{matrix}\right] modelx=⎣⎡0000−10001⎦⎤modely=⎣⎡0000010−10⎦⎤

得到 Z x , Z y Z_x\ ,Z_y Zx,Zy；

接着根据公式对图像的所有点进行处理： Z = Z x 2 + Z y 2 Z=\sqrt{Z_x^2+Z_y^2} Z=Zx2+Zy2 ；

2） S o b e l Sobel Sobel算子：先做加权平滑处理，再做微分运算，平滑部分的权值有些差异，对噪声具有一定的抑制能力，但不能完全排除虚假边缘。虽然这两个算子边缘定位效果不错，但检测出的边缘容易出现多像素宽度。

首先使用以下模板分别对图像进行卷积操作：

m o d e l x = [ − 1 − 2 − 1 0 0 0 1 2 1 ] m o d e l y = [ 1 0 − 1 2 0 − 2 1 0 − 1 ] modelx=\left[\begin{matrix}-1&-2&-1\\0&0&0\\1&2&1\\\end{matrix}\right]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ modely=\left[\begin{matrix}1&0&-1\\2&0&-2\\1&0&-1\\\end{matrix}\right] modelx=⎣⎡−101−202−101⎦⎤modely=⎣⎡121000−1−2−1⎦⎤

得到 Z x , Z y Z_x\ ,Z_y Zx,Zy；

接着根据公式对图像的所有点进行处理： Z = Z x 2 + Z y 2 Z=\sqrt{Z_x^2+Z_y^2} Z=Zx2+Zy2 ；

3） P r e w i t t Prewitt Prewitt算子

首先使用以下模板分别对图像进行卷积操作：

m o d e l x = [ − 1 − 1 − 1 0 0 0 1 1 1 ] m o d e l y = [ 1 0 − 1 1 0 − 1 1 0 − 1 ] modelx=\left[\begin{matrix}-1&-1&-1\\0&0&0\\1&1&1\\\end{matrix}\right]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ modely=\left[\begin{matrix}1&0&-1\\1&0&-1\\1&0&-1\\\end{matrix}\right] modelx=⎣⎡−101−101−101⎦⎤modely=⎣⎡111000−1−1−1⎦⎤

得到 Z x Z y Z_x \,Z_y ZxZy；

接着根据公式对图像的所有点进行处理： Z = Z x 2 + Z y 2 Z=\sqrt{Z_x^2+Z_y^2} Z=Zx2+Zy2 ；

4） L O G LOG LOG算子：先用高斯函数平滑滤波，后用 L a p l a c i a n Laplacian Laplacian算子检测边缘，克服 L a p l a c i a n Laplacian Laplacian算子抗噪声差的缺点，但同时也平滑掉比较尖锐的边缘，尖锐边缘无法检被测到。

使用以下模板对图像进行卷积操作：

m o d e l = [ − 2 − 4 − 4 − 4 − 2 − 4 0 8 0 − 4 − 4 8 24 8 − 4 − 4 0 8 0 − 4 − 2 − 4 − 4 − 4 − 2 ] model=\left[\begin{matrix}-2&-4&-4&-4&-2\\-4&0&8&0&-4\\-4&8&24&8&-4\\-4&0&8&0&-4\\-2&-4&-4&-4&-2\\\end{matrix}\right] model=⎣⎢⎢⎢⎢⎡−2−4−4−4−2−4080−4−48248−4−4080−4−2−4−4−4−2⎦⎥⎥⎥⎥⎤

通过以上四种算子对图像进行处理可以得到图像的边缘信息，而锐化图像 = 原图像 + 加重的边缘，由此可得到根据四种算子处理的锐化图像。

7. 实验中直接使用的函数

（1）imread(path)函数

从图像所在路径读取图像的数据信息存为矩阵。

（2）imshow(image)函数

将读取到的图像显示到figure中。

（3）subplot(m,n,p)函数

subplot函数是将多个图画到一个平面上的工具。其中m、n表示一个m行n列的大画框，可显示 m ∗ n m*n m∗n个图，p表示图所在位置。

（4）构造函数

function[输出形参]=函数名([输入形参])

函数体

（5）size()函数

[m,n] = size(X)

返回矩阵X的尺寸信息，并存储在m、n中。其中m中存储的是行数，n中存储的是列数。

（6）zeros(m,n)函数

产生 m ∗ n m*n m∗n的double类型零矩阵。

（7）double(x)

增加数据的精度，使计算更加准确。

（8）imnoise(I,type)

向亮度图I中添加指定类型的噪声。type是字符串，可以是以下值。

‘‘gaussian’’(高斯噪声);’‘localvar’’(均值为零，且一个变量与图像亮度有关);’‘poisson’’(泊松噪声)；’‘salt&pepper’’(椒盐噪声)；’‘speckle’’(乘性噪声)

（9）sum(A):求和

（10）sort(A):排序

（11）sqrt(x):开根号

（12）strcmp(choice,’’):比较字符串是否相等

3.2 实验代码

【sy2.m】

clear all;close all;clc;%% 1、为图像分别加高斯噪声和椒盐噪声，采用中值滤波方法对受噪声干扰的图像滤波，窗口分别采用3*3,5*5,7*7。I=imread('Fig1006(a).tif');%读取灰度图像J=imnoise(I,'salt & pepper',0.06); %给灰度图像添加椒盐噪声K=imnoise(I,'gaussian',0.05,0.05); %给灰度图像添加高斯噪声 figure,suptitle('原图+加噪声的图像'); %Figure 1subplot(131), imshow(I), title('原图'); subplot(132), imshow(J), title('椒盐噪声');subplot(133), imshow(K), title('高斯噪声');J1=mymiddle(J,3);J2=mymiddle(J,5);J3=mymiddle(J,7);figure,suptitle('中值滤波处理椒盐噪声'); %Figure 2subplot(131), imshow(J1), title('3×3 模板的中值滤波');subplot(132), imshow(J2), title('5×5 模板的中值滤波');subplot(133), imshow(J3), title('7×7 模板的中值滤波');K1=mymiddle(K,3);K2=mymiddle(K,5);K3=mymiddle(K,7);figure,suptitle('中值滤波处理高斯噪声');%Figure 3subplot(131), imshow(K1), title('3×3 模板的中值滤波');subplot(132), imshow(K2), title('5×5 模板的中值滤波');subplot(133), imshow(K3), title('7×7 模板的中值滤波');%% 2、对受噪声干扰的图像进行均值滤波。J4=myaverage(J,3);K4=myaverage(K,3);figure,suptitle('3×3 模板的均值滤波'); %Figure 4subplot(121), imshow(J4), title('均值滤波处理椒盐噪声');subplot(122), imshow(K4), title('均值滤波处理高斯噪声');%% 3、分别采用Roberts 算子、Sobel算子、Prewitt算子和Log算子进行图像锐化。I=imread('Fig1006(a).tif');S1=mysharpen(I,'Roberts');S2=mysharpen(I,'Sobel');S3=mysharpen(I,'Prewitt');S4=mysharpen(I,'Log');figure,suptitle('图像锐化[1]'); %Figure 5subplot(131), imshow(I), title('原图');subplot(132), imshow(I+S1), title('Roberts算子');subplot(133), imshow(I+S2), title('Sobel算子');figure,suptitle('图像锐化[2]'); %Figure 6subplot(121), imshow(I+S3), title('Prewitt算子');subplot(122), imshow(I+S4), title('Log算子');

【mymiddle.m】

% 函数mymiddle:根据模板对图像进行中值滤波% 输入参数:image:原图像%k:模板大小为k×k（k为奇数）% 输出参数:中值滤波后的图像数据OUT% 使用函数:size(x):求矩阵大小% double(x):增加精度% zeros():为矩阵分配空间% sort(A):排序function OUT=mymiddle(image,k)HW=size(image);%获取原图像的大小image=double(image);%运算的时候用double,为了保持运算精度wall=(k-1)/2;nHW=[HW(1)+2*wall,HW(2)+2*wall];%加上墙（补充的边沿）之后的数据区大小data=zeros(nHW);%为数据区分配空间data(wall+1:wall+HW(1),wall+1:wall+HW(2))=image(:,:);%为数据区赋值%加'墙'for i = 1:walldata(wall+1:wall+HW(1),i)=data(wall+1:wall+HW(1),wall+1);%左data(wall+1:wall+HW(1),i+(nHW(2)-wall))=data(wall+1:wall+HW(1),wall+HW(2));%右endfor i = 1:walldata(i,:)=data(wall+1,:);%上data(i+(nHW(1)-wall),:)=data(nHW(1)-wall,:);%下end%取中值for i = wall+1:wall+HW(1)for j = wall+1:wall+HW(2)A=zeros(1,k*k);%为k×k模板的数据分配空间for m = i-wall:i+wallfor n = j-wall:j+wallA((m-i+wall)*k+(n-j+wall+1))=data(m,n);endendA=sort(A);data(i,j)=A((k*k+1)/2);endendOUT=uint8(data(wall+1:wall+HW(1),wall+1:wall+HW(2)));%存储的时候的一般存储为uint8类型，节省存储空间,显示的时候用uint8

【myaverage.m】

% 函数myaverage:根据模板对图像进行中值滤波% 输入参数:image:原图像%k:模板大小为k×k（k为奇数）% 输出参数:均值滤波后的图像数据OUT% 使用函数:size(x):求矩阵大小% double(x):增加精度% zeros():为矩阵分配空间% sum(A):求和function OUT=myaverage(image,k)HW=size(image);%获取原图像的大小image=double(image);%运算的时候用double,为了保持运算精度wall=(k-1)/2;nHW=[HW(1)+2*wall,HW(2)+2*wall];%加上墙（补充的边沿）之后的数据区大小data=zeros(nHW);%为数据区分配空间data(wall+1:wall+HW(1),wall+1:wall+HW(2))=image(:,:);%为数据区赋值%加'墙'for i = 1:walldata(wall+1:wall+HW(1),i)=data(wall+1:wall+HW(1),wall+1);%左data(wall+1:wall+HW(1),i+(nHW(2)-wall))=data(wall+1:wall+HW(1),wall+HW(2));%右endfor i = 1:walldata(i,:)=data(wall+1,:);%上data(i+(nHW(1)-wall),:)=data(nHW(1)-wall,:);%下end%求均值for i = wall+1:wall+HW(1)for j = wall+1:wall+HW(2)A=zeros(1,k*k);%为k×k模板的数据分配空间for m = i-wall:i+wallfor n = j-wall:j+wallA((m-i+wall)*k+(n-j+wall+1))=data(m,n);endenddata(i,j)=sum(A)/(k*k);%求和/模板大小=均值endendOUT=uint8(data(wall+1:wall+HW(1),wall+1:wall+HW(2)));

【mysharpen.m】

% 函数mysharpen:根据选择的算子对图像进行锐化处理% 输入参数:I:原图像%modelx:水平边缘算子%modely:垂直边缘算子% 输出参数:锐化处理后的图像OUT% 使用函数:size(x):求矩阵大小% sqrt(x):开根号% strcmp(choice,''):比较字符串是否相等% myconvolve(I，model):卷积操作function OUT=mysharpen(I,choice)[x,y]=size(I);%% 采用Roberts算子进行边缘检测if strcmp(choice,'Roberts')modelx=[0,0,0;0,-1,0;0,0,1];%水平边缘算子modely=[0,0,0;0,0,-1;0,1,0];%垂直边缘算子b=myconvolve(I,modelx);c=myconvolve(I,modely);for i=1:xfor j=1:yb(i,j)=sqrt(b(i,j)^2+c(i,j)^2); endendOUT=uint8(b);end%% 采用Sobel算子进行边缘检测if strcmp(choice,'Sobel')modelx=[-1,-2,-1;0,0,0;1,2,1];%水平边缘算子modely=[1,0,-1;2,0,-2;1,0,-1];%垂直边缘算子b=myconvolve(I,modelx);c=myconvolve(I,modely);for i=1:xfor j=1:yb(i,j)=sqrt(b(i,j)^2+c(i,j)^2); endendOUT=uint8(b);end%% 采用Prewitt算子进行边缘检测if strcmp(choice,'Prewitt')modelx=[-1,-1,-1;0,0,0;1,1,1];%水平边缘算子modely=[-1,0,1;-1,0,1;-1,0,1];%垂直边缘算子b=myconvolve(I,modelx);c=myconvolve(I,modely);for i=1:xfor j=1:yb(i,j)=sqrt(b(i,j)^2+c(i,j)^2); endendOUT=uint8(b);end%% 采用Log算子进行边缘检测if strcmp(choice,'Log')model=[-2,-4,-4,-4,-2;-4,0,8,0,-4;-4,8,24,8,-4;-4,0,8,0,-4;-2,-4,-4,-4,-2];%模板b=myconvolve(I,model);OUT=uint8(b);end

【myconvolve.m】

% 函数myconvolve:根据模板对图像进行卷积操作% 输入参数:I:原图像%modelx:模板% 输出参数:卷积操作后的图像数据OUT% 使用函数:size(x):求矩阵大小% double(x):增加精度% zeros():为矩阵分配空间function OUT=myconvolve(image,model)HW=size(image); %获取原图像的大小image=double(image);[x,y]=size(model);%模板大小OUT=zeros(HW); %为输出分配空间%% 加边沿wall=(x-1)/2;nHW=[HW(1)+2*wall,HW(2)+2*wall];%加上墙（补充的边沿）之后的数据区大小data=zeros(nHW);%为数据区分配空间data(wall+1:wall+HW(1),wall+1:wall+HW(2))=image(:,:);%为数据区赋值for i = 1:walldata(wall+1:wall+HW(1),i)=data(wall+1:wall+HW(1),wall+1);%左data(wall+1:wall+HW(1),i+(nHW(2)-wall))=data(wall+1:wall+HW(1),wall+HW(2));%右endfor i = 1:walldata(i,:)=data(wall+1,:);%上data(i+(nHW(1)-wall),:)=data(nHW(1)-wall,:);%下end%% 计算for i = wall+1:wall+HW(1)for j = wall+1:wall+HW(2)for m = 1:xfor n = 1:yOUT(i-wall,j-wall) = OUT(i-wall,j-wall) + data(i-wall+m-1,j-wall+n-1)*model(m,n);%对整幅图像进行卷积操作endendendend