Spatial Filtering and Median Filtering

图像的代数运算是图像的标准算术操作的实现方法，是两幅输入图像之间进行点对点的加、减、乘、除运算后得到输出图像的过程。空域滤波是在图像空间中借助模板对图像进行领域操作,在待处理图像中逐点地移动掩膜,并在掩膜扫过的区域进行指定运算，其理论基础是空间卷积，目的是改善影像质量。平滑滤波用于模糊处理、减小噪声，图像中每个像素用3X3区域内像素均值代替，以减小图像灰度的尖锐变化。拉普拉斯是一种微分算子，它的应用强调图像中灰度的突变及降低灰度慢变化的区域。将原始图像和拉普拉斯图像叠加加在一起的简单方法可以保护拉普拉斯锐化处理效果，同时又能复原背景信息。中值滤波的响应基于图像滤波器包围的图像区域中像素的排序，然后用排序结果中间值代替中心像素的值，此种方法对于一定类型的随机噪声相当有效，如椒盐噪声。本实验目的注重于图像的加、减、乘、除运算以及各种空域滤波的实现

一、技术论述
1.算术运算
（1）图像加法运算
图像加法运算一般用于对同一场景的多幅图像求平均效果，以便有效地降低具有叠加性质的随机噪声，于那些经过长距离模拟通讯方式传送的图像，这种处理是必不可少的。两幅图像f(x,y)与h(x,y)的和可以表示为：g(x,y)=f(x,y)+h(x,y)
（2）图像减法运算
图像减法也称为差分方法，是一种常用于检测图像变化及运动物体的图像处理方法。可以用来检测一系列相同场景图像的差异。图像的差异是通过计算这两幅图像所有对应像素点的差而得出的。减法处理最主要的作用就是增强两幅图像的差异两幅图像。f(x,y)与h(x,y)的差异表示为：g(x,y)=f(x,y)-h(x,y)
（3）图像乘法运算
两幅图像进行乘法运算可以实现掩膜操作，即屏蔽掉图像的某些部分。一幅图像乘以一个常数通常被称为缩放，如果使用的缩放因数大于1，那么将增强图像的亮度；如果缩放因数小于1，则会使图像变暗。两幅图像f(x,y)与h(x,y)的积可以表示为：g(x,y)=f(x,y)*h(x,y)、g(x,y)=k*f(x,y)
（4）图像除法运算
除法运算可用于校正成像设备的非线性影响，也可以用来检测两幅图像的区别。两幅图像相除看做一幅的取反图像与另一幅图像相乘。除法操作给出的相应像素值的变化比率，而不是每个像素的绝对差异，因而图像除法也称为比率变化。两幅图像f(x,y)与h(x,y)相除可以表示为：g(x,y)=f(x,y)/h(x,y)
2.空域滤波
空间滤波处理就是再待处理图像中逐点地移动掩膜。在每一点（x,y）处，滤波器在该店的响应通过事先定义的关系计算。对于线性空间滤波，其响应由滤波器系数与滤波掩膜扫过区域的相应像素的乘积之和给出。如图所示为3*3的掩膜，在图像中的点(x,y)处，用该掩膜线性滤波的响应R为：
R = w(-1,-1)f(x-1,y-1)+ w(-1,0)f(x-1,y) +…+ w(0,0)f(x,y) +
+ w(1,0)f(x+1,y) + w(1,1)f(x+1,y+1)

一般3*3空间滤波掩膜
系数w(0,0)与图像值f(x,y)相符合是指当乘积求和计算发生时，掩膜位于(x,y)的中心。
空域滤波基本思路：
a) 设计模版类型、模板参数和模板运算方法；
b) 将模板中心与图像中某一像素位置重合；
c) 对模板覆盖的图像区域中的像素点进行运算；
d) 将运算结果替代图像中模板中心对应位置的像素；

（1）平滑滤波
平滑滤波器用于模糊处理和降低噪声。模糊处理经常用于预处理任务中，例如在（大）目标提取之前去除图像中的一些琐碎细节，以及其桥接直线或曲线的缝隙。通过线性滤波和非线性滤波模糊处理，可以降低噪声。
平滑滤波一般有两种掩膜，直接平均掩膜与加权平均掩膜，其结构对应下图：

两个3X3均值滤波器掩膜

（2）中值滤波
中值滤波是将邻域内像素（在中值计算中包括原像素的值）灰度的中值代替该像素的值，中值滤波应用非常普遍，，其对于一定类型的随机噪声有很好的去噪能力，比小尺度的线性平滑滤波器的模糊程度明显要低。尤其是处理椒盐噪声有非常好的效果。
中值滤波的思路是针对某个像素f(x,y)，取其3X3邻域内的9个像素，将这九个像素点按从小到大的顺序排列，并取中间值作为像素点f(x,y)的灰度值。

（3）二阶微分图像增强——拉普拉斯算子
拉普拉斯是一种微分算子，应用强调图像中灰度的突变及降低灰度慢变化的区域，这将产生衣服把图像中的浅灰色边线和突变点叠加到暗背景中的图像。将原始图像和拉普拉斯图像叠加加在一起的简单方法可以保护拉普拉斯锐化处理效果，同时又能复原背景信息。
拉普拉斯定义为下式：

一般二维拉普拉斯数字实现见下式：

实验中使用的两种改进拉普拉斯算子掩膜结构

所以拉普拉斯变换对图像增强的基本方法可以表示为下式：

当拉普拉斯掩膜中心系数为负时选择相减，为正时选择相加。

（4）一阶微分图像增强——梯度法
邻域平均法或加权平均法可以平滑图像，反过来利用对应的微分算法可以锐化图像。微分算法是求信号的变化率，有加强高频分量的作用，从而使图像轮廓清晰。由于图像模糊胡实质是图像受到平均或积分运算造成的，所以为了把图像中任何方向伸展的边缘的模糊的轮廓变得清晰，可以对图像进行逆运算如微分运算，从而使图像清晰化。在图像处理中，一阶微分是通过梯度算法来实现的。对于一幅图像用函数f(x,y)表示，定义在f(x,y)在点(x,y)处的梯度是一个矢量，定义为：

梯度的方向在函数f(x,y)最大变化率的方向上，梯度的膜值G[f(x,y)]可以由以下公式算出：

对于数字图像而言，微分可用差分来近似。因此上式可以写成：

于此，本实验采用如下图所示的梯度法的掩膜，其结构如下：

则在点w5使用上述掩膜，其拉普拉斯变换结果近似等于：
| (w7+2w8+w9）-(w1+2w2+w3) |+ | (w3+2w6+w9)-(w1+2w4+w7)|

二、实验结果讨论
1.图像的代数运算
由于图像的灰度级处于[0,255]，所以超出这个范围后像素值会被置0或255。
(1)图像相加一般用于多幅图像求平均效果，以便有效地降低具有叠加性质的随机噪声，见图1.1。
(2)图像相减可以显示出两幅图之间的差异，见图1.1。如果需要可以将结果加到原图上，可以突出差异细节。
(3)图像相乘实现掩模操作，即屏蔽掉图像的某些部分。一幅图像乘以一个常数被称为缩放，缩放将产生比简单添加像素偏移量自然得多的明暗效果，因为这种操作能够更好的维持图像的相关对比度。见图1.2。
(4)图像相除用于校正成像设备的非线性影响，也可以用来检测两幅图像间的区别，但是除法操作给出的是相应像素值的变化比率，而不是每个像素的绝对差异，因而图像除法也成为比率变换。见图1.2，只有差异特别大的地方被显示出来，如想看到更多细节，也可在相除后乘上一个因数提高灰度差异不是特别大的地方的灰度值。

2.平滑滤波
平滑滤波用于模糊处理、减小噪声，图像中每个像素用3X3区域内像素均值代替，减小了图像灰度的尖锐变化区域。见图2，相对于原图，处理后图像边缘明显模糊。

3.拉普拉斯算子滤波法
微分运算能够强调图像中灰度的突变及降低灰度慢变化的区域，见图3.1，原图中的边缘部分明显被增强，图3.2是另外一种系数的掩膜滤波后的结果，边缘突出更加明显。

4.梯度法滤波
梯度法属于一阶微分法，能加强高频分量的作用，从而使图像轮廓清晰，见图4.1，原图的轮廓被提取了出来，再加到原图便可以起到增强图片边缘的作用，见图4.2。

5.中值滤波
中值滤波对于去除一些随机噪声有非常好的效果，如椒盐噪声，中值滤波使拥有不同灰度的点看起来更接近它的邻近值，去除那些邻域中相对较暗、亮的孤立点。图5.1为添加椒盐噪声后的图像，图5.2为中值滤波后的图像。对于课本图3.37(b)所示的平滑滤波结果图，虽然可见噪声减少了，但是付出的代价是图像变得模糊了。

三、实验结果

图1.1

图1.2

图2

图3.1

图3.2

图4.1

图4.2

图5.1

图5.2

附录
1.图像加、减、乘、除程序：
function daishuyunsuan()
i=imread('Fig1.10(4).jpg');
j=imread('Fig1.10(5).jpg');
subplot(2,2,1),imshow(i); title('原始图像i')
subplot(2,2,2),imshow(j); title('原始图像j')
g=imsubtract(i,j);
subplot(2,2,3),imshow(uint8(g));title('图像i-k')
g=imadd(i,j);
subplot(2,2,4),imshow(uint8(g)); title('图像i+k')
figure
g=immultiply(i,j);
subplot(2,2,1),imshow(uint8(g));title('图像i*k')
g=immultiply(i,2);
subplot(2,2,2),imshow(uint8(g));title('图像i*常数')
g=imdivide(i,j);
subplot(2,2,3),imshow(uint8(g));title('图像i/k')
g=imdivide(i,2);
subplot(2,2,4),imshow(uint8(g));title('图像i/常数')

2.平滑滤波程序
function pinghua(w)
%w=[1 1 1;1 1 1; 1 1 1]*1/9;
im=imread('Fig3.38(a).jpg');
figure
subplot(1,2,1);
imshow(im);
im=double(im);
[x,y]=size(im)
% X=im;
for i=2:(x-1)
    for j=2:(y-1)
     im1(i,j)=im(i-1,j-1)*w(1,1)+im(i-1,j)*w(1,2)+im(i-1,j+1)*w(1,3)+im(i,j-1)*w(2,1)+im(i,j)*w(2,2)+im(i,j+1)*w(2,3)+im(i+1,j-1)*w(3,1)+im(i+1,j)*w(3,2)+im(i+1,j+1)*w(3,3);%+im(i,j);   
    end
end

[xx,yy]=size(im1)
subplot(1,2,2);
imshow(uint8(im1))



3.拉普拉斯算子滤波程序
function lapulasi(w)
%w=[-1 -1 -1;-1 8 -1; -1 -1 -1];
im=imread('Fig3.38(a).jpg');
figure
% title('拉普拉斯算子图像增强')
subplot(1,2,1);
imshow(im);
% title('w=[-1 -1 -1;-1 8 -1; -1 -1 -1];
im=imread('Fig3.38(a).jpg');
figure
% title('拉普拉斯算子图像增强')
subplot(1,2,1);
imshow(im);
% title('原图像')
im=double(im);
[x,y]=size(im)
% X=im;
for i=2:(x-1)
    for j=2:(y-1)
   im1(i,j)=im(i-1,j-1)*w(1,1)+im(i-1,j)*w(1,2)+im(i-1,j+1)*w(1,3)+im(i,j-1)*w(2,1)+im(i,j)*w(2,2)+im(i,j+1)*w(2,3)+im(i+1,j-1)*w(3,1)+im(i+1,j)*w(3,2)+im(i+1,j+1)*w(3,3)+im(i,j);   
    end
end
[xx,yy]=size(im1)
subplot(1,2,2);
imshow(uint8(im1))')

4.梯度法滤波程序
function  tidu(w,u)
%w=[-1 -2 -1;0 0 0; 1 2 1];
%u=[-1 0 1;-2 0 2;-1 0 1];
im=imread('Fig3.38(a).jpg');
figure
% title('拉普拉斯算子图像增强')
subplot(1,2,1);
imshow(im);
% title('原图像')
im=double(im);
[x,y]=size(im)
% X=im;
for i=2:(x-1)
    for j=2:(y-1)
im1(i,j)=abs(im(i-1,j-1)*w(1,1)+im(i-1,j)*w(1,2)+im(i-1,j+1)*w(1,3)+im(i,j-1)*w(2,1)+im(i,j)*w(2,2)+im(i,j+1)*w(2,3)+im(i+1,j-1)*w(3,1)+im(i+1,j)*w(3,2)+im(i+1,j+1)*w(3,3))   +abs(im(i-1,j-1)*u(1,1)+im(i-1,j)*u(1,2)+im(i-1,j+1)*u(1,3)+im(i,j-1)*u(2,1)+im(i,j)*u(2,2)+im(i,j+1)*u(2,3)+im(i+1,j-1)*u(3,1)+im(i+1,j)*u(3,2)+im(i+1,j+1)*u(3,3))+im(i,j);   
    end
end

[xx,yy]=size(im1)
subplot(1,2,2);
imshow(uint8(im1))

5.添加椒盐噪声程序
function noise()
im=imread('Fig5.07(a).jpg');
figure
im1=imnoise(im,'salt & pepper',0.2);
subplot(1,2,1);
imshow(im);
subplot(1,2,2);
imshow(uint8(im1))

6.中值滤波程序
function zhongzhilvbo()
im=imread('jioayanjieguo.jpg');
figure
subplot(1,2,1);
imshow(im);
im=double(im);
[x,y]=size(im)
% X=im;
for i=2:(x-1)
    for j=2:(y-1)
        w=[im(i-1,j-1) im(i-1,j) im(i-1,j+1) im(i,j-1) im(i,j) im(i,j+1) im(i+1,j-1) im(i+1,j) im(i+1,j+1)];
        w=sort(w);
        im1(i,j)=w(5);   
    end
end

[xx,yy]=size(im1)
subplot(1,2,2);
imshow(uint8(im1))