详解数字图像的滤波和边缘检测

一、图像滤波

图像滤波的主要目的就是在尽量保留图像细节特征的条件下对目标图像的噪声进行抑制。图像滤波其主要分为线性滤波器和非线性滤波器。

1.1 线性滤波器

1.1.1 均值滤波

原理：在图像上，对待处理的像素给定一个模板，该模板包括了其周围的邻近像素。将模板中的全体像素的均值来替代原来的像素值的方法。以 \(3\times3\) 的均值滤波器的模板为：

\[w= \begin{bmatrix}
1&1&1\\
1&1&1\\
1&1&1
\end{bmatrix}
\]

模板所覆盖的图像的像素为：

\[Image = \begin{bmatrix}
(i-1,j-1)&(i-1,j)&(i-1,j+1)\\
(i,j-1)&(i,j)&(i,j+1)\\
(i+1,j-1)&(i+1,j-1)&(i+1,j+1)
\end{bmatrix}
\]

因此，模板中心点所对应的像素值可以表示为：

\[g(i,j) = \dfrac{1}{N} \sum_{i=1}^{N}Image \odot w = \dfrac{1}{N}\sum_{i=1}^{N}\sum_{m=-1}^{1}\sum_{n=-1}^{1}f(i+m,j+n)
\]

在进行滤波之前，一般会先对原图像进行填充，填充的像素大小为：\(r = floor(模板宽度/2)\)，不然的话始终无法对图像边缘的像素进行滤波。

均值滤波的特点：均值滤波能够降低图像中的尖锐变化，去除图像中的噪声，且滤波窗口越大，去噪效果越好同时图像也会变得越模糊

1.1.1 高斯滤波

原理：其实高斯滤波和均值滤波的原理差不多，都是用一个模板（或称卷积、掩模）扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值。区别就在于，高斯滤波中模板的数值是通过高斯函数生成的。还是以 \(3\times3\) 的模板为例，以模板的左上角为坐标原点，模板中心的坐标为 \((1,1)\)，则生成的高斯模板为：

\[w(i,j) = \dfrac{1}{2\pi \sigma^{2}}exp{-\dfrac{(i-1)^2+(j-1)^2}{2\sigma^{2}}}
\]

因此，模板中心点所对应的像素值可以表示为：

\[g(i,j) = Image \odot \begin{bmatrix}
w_{(i-1,j-1)}&w_{(i-1,j)}&w_{(i-1,j+1)}\\
w_{(i,j-1)}&w_{(i,j)}&w_{(i,j+1)}\\
w_{(i+1,j-1)}&w_{(i+1,j)}&w_{(i+1,j+1)}
\end{bmatrix} = \sum_{m=-1}^{1}\sum_{n=-1}^{1} w(i+m,j+n)f(i+m,j+n)
\]

高斯滤波的特点：根据公式可以得到，在计算高斯滤波模板时，高斯滤波模板的值由中间向四周递减，且标准差越小，二维高斯图像越窄小，平滑效果不明显；标准差越大，而为高斯图像越矮宽，滤波效果比较明显。下图展示了不同 \(\sigma\) 的\(7\times7\) 高斯模板的差别：

1.2 非线性滤波器

1.2.1 中值滤波器

原理：中值滤波器就是取模板覆盖区域的排序之后的中间值作为该模板区域内中心的像素值。还是以 \(3\times 3\) 中值滤波器模板为例，此时覆盖的像素为：

\[Image = \begin{bmatrix}
(i-1,j-1)&(i-1,j)&(i-1,j+1)\\
(i,j-1)&(i,j)&(i,j+1)\\
(i+1,j-1)&(i+1,j-1)&(i+1,j+1)
\end{bmatrix}
\]

那么，中值滤波后的像素值为：

\[g(i,j) = median(Image)
\]

中值滤波器的特点：中值滤波器很适合处理椒盐噪声。

1.2.2 双边滤波器

设计原因：在均值滤波器和高斯滤波器时，滤波器的模板值仅仅考虑到了像素的空间位置，潜在的假设为：离模板中心点越近，那么他们的像素值也越接近。但是在纹理比较复杂的区域，这种假设很难成立，因此采用上述滤波器进行滤波时，图像的边缘信息保持不够理想。那么能不能设计一款滤波器，能够识别纹理复杂的区域，从而实现更好的边缘保持i效果？

原理：双边滤波器分为两个部分，首先采用高斯滤波器生成只跟像素位置相关的模板，然后根据像素值，生成另外一个跟像素值相关的模板。两个模板逐像素点击就形成了双边滤波模板。具体描述如下：

\[w1(i,j,k) = exp(-\dfrac{(i-k)^2+(j-k)^2}{2\sigma_d^{2}}) \\
w2(i,j,k) = exp(-\dfrac{||f(i,j)-f(k,k)||^{2}}{2\sigma_d^{2}})
\]

因此，最终的双边滤波模板为：

\[w = w1 * w2 = exp(-\dfrac{(i-k)^2+(j-k)^2}{2\sigma_d^{2}}-\dfrac{||f(i,j)-f(k,k)||^{2}}{2\sigma_d^{2}})
\]

在双边滤波器中多考虑了值域模板，当两个位置像素越接近，那么 \(w2\) 对应位置的权重就越大，否则就会越小。因此通过 \(w2\) 我们在滤波器中保护了图像的边缘信息。其模板的图像如下：

双边滤波器的特点：双边滤波器是一种可以保边缘信息的去噪滤波器。

1.3 滤波器的 Matlab 代码实现

clear all
clc
image = imread(\'../测试图像/2.bmp\');
[H,W,C] = size(image);  %记录读入图像的长宽

H_win = 7;  %定义模板窗口的宽度
image_input = double(rgb2ycbcr(image)); %转换成Ycbcr颜色空间
image_input = image_input(:,:,1); %取Y通道来进行处理

paddle = floor(H_win/2);    %计算需要填充的像素宽度
image_data = zeros(H+2*paddle, W+2*paddle); %0填充

Start = paddle+1;   %计算图像的开始和结束位置
End = H+paddle;
image_data(Start:End, Start:End) = image_input(:,:);%将图像放在填充图像的中间

%生成高斯滤波模板
sigma_g = 3;
sigma_r = 6;
[x,y] = meshgrid(-paddle:paddle,-paddle:paddle);
k_gaussian1=exp(-(x.^2+y.^2)/(2*sigma_g^2));     %以距离作为自变量高斯滤波器

k_gaussian=k_gaussian1/sum(k_gaussian1(:));%归一化
Map_X = 1:H_win; Map_Y = 1:H_win; %作3维图的横坐标
surf(Map_X,Map_Y,k_gaussian)    %绘制高斯滤波核的图


for i=Start:1:End
    for j=Start:1:End
        image_window = image_data(i-paddle:i+paddle, j-paddle:j+paddle);
        
        Mean_image(i-paddle,j-paddle) = mean(mean(image_window));   %均值滤波
        Median_image(i-paddle,j-paddle) = median(image_window(:));  %中值滤波
        
        Gaussian_image(i-paddle, j-paddle) = sum(sum(image_window.*k_gaussian));  %高斯滤波
        
        double_w = exp(-(image_window-image_data(i,j)).^2/(2*sigma_r^2)); %以周围和当前像素灰度差值作为自变量的高斯滤波器
        double_w = k_gaussian1 .* double_w;
        double_w = double_w/sum(double_w(:));
        
        Double_image(i-paddle, j-paddle) = sum(sum(image_data(i-paddle:i+paddle, j-paddle:j+paddle).*double_w));  %双边滤波
    end
end
figure, imshow(uint8(image_input)), title(\'原图\')
subplot(2,2,1), imshow(uint8(Mean_image)), title(\'5*5均值\')
subplot(2,2,2), imshow(uint8(Median_image)), title(\'中值滤波\')
subplot(2,2,3), imshow(uint8(Gaussian_image)), title(\'5*5,1.5 高斯\')
subplot(2,2,4), imshow(uint8(Double_image)), title(\'5*5,1.5 双边\')

代码输出结果如下：

二、图像边缘检测

边缘检测是为了将其周围像素灰度有阶跃变化的像素检测出来，这些像素组成的集合就是该图像的边缘。比较常用的边缘检测方法就是考察每个像素在某个领域内灰度的变化，然后利用边缘临近一阶或二阶方向导数变化规律检测边缘，即边缘检测局部算法

2.1 一阶边缘检测算子

2.1.1 Sobel 算子

原理：Sobel 算子包括两组 \(3\times3\) 的矩阵，左边的表示垂直，右边的表示水平。将它与图像作平面卷积，即可分别得出垂直及水平的亮度差分近似值。

\[G_x = \begin{bmatrix} -1&0&1\\-2&0&2\\-1&0&1\\\end{bmatrix}, G_y = \begin{bmatrix} 1&2&1\\0&0&0\\-1&-2&-1\\\end{bmatrix}
\]

假设 \(O_x\) 和 \(O_{y}\) 分别用来表示被模板 \(G_x\) 和 \(G_{y}\) 卷积后的结果，那么最终的输出图像边缘可以表示为：

\[O = \sqrt{O_x^{2}+O_y^{2}}
\]

Sobel算子的特点：Soble算子在水平和垂直两个方向上求导，得到的是图像在X方法与Y方向梯度图像。但是比较敏感，容易受影响，要通过高斯模糊（平滑）来降噪。

2.1.2 Canny 算子

原理：Canny算子是由计算机科学家John F. Canny于1986年提出的一种边缘检测算子，是目前理论上相对最完善的一种边缘检测算法。主要包括以下几个步骤：

高斯滤波对图像进行预处理，提前弱化噪声，增加边缘提取准确度
Sobel 算子进行像素梯度计算，得到梯度强度矩阵
非极大值抑制：非极大值像素梯度抑制起到将边缘“瘦身”的作用。其基本方法是将当前像素梯度强度与沿正负梯度方向上的相邻像素的梯度强度进行比较，若其最大（即为极值），则保留该像素为边缘点，若不是最大，则对其进行抑制，不将其作为边缘点。
滞后阈值处理：定义一个高阈值和一个低阈值。梯度强度低于低阈值的像素点被抑制，不作为边缘点；高于高阈值的像素点被定义为强边缘，保留为边缘点；处于高低阈值之间的定义为弱边缘，留待进一步处理。
孤立弱边缘抑制：通常而言，由真实边缘引起的弱边缘像素点将连接到强边缘像素点，而噪声响应则未连接。通过查看弱边缘像素及其8个邻域像素，可根据其与强边缘的连接情况来进行判断。一般，可定义只要其中邻域像素其中一个为强边缘像素点，则该弱边缘就可以保留为强边缘，即真实边缘点。

Canny检测的特点：Canny 算子对于滤波参数和高、低阈值的选取还是较为敏感，使得实际应用过程中需要反复调试。

2.2 二阶边缘检测算子

2.2.1 拉普拉斯算子

原理：拉普拉斯算子是二阶微分算子，实际上就是图像梯度的散度：

\[\nabla^{2}f = \dfrac{\partial^2f}{\partial x^2}+\dfrac{\partial^2f}{\partial y^2}
\]

然而：

\[\dfrac{\partial^2f}{\partial x^2}=f(x,y+1)-2f(x,y)+f(x,y-1) \\
\dfrac{\partial^2f}{\partial y^2}=f(x+1,y)-2f(x,y)+f(x-1,y)
\]

所以，卷积的模板可以写为：

\[\nabla^{2}f = \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}, \nabla^{2}f = \begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix}
\]

拉普拉斯算子的特点：Laplacian算子法对噪声比较敏感，所以很少用该算子检测边缘，而是用来判断边缘像素视为与图像的明区还是暗区。Laplacian算子是各向同性的，能对任何走向的界线和线条进行锐化，无方向性。

2.3 边缘检测的 Matlab 代码实现

clear all
clc
image = imread(\'../测试图像/2.bmp\');
[H,W,C] = size(image);  %记录读入图像的长宽

H_win = 21;  %定义模板窗口的宽度
image_input = double(rgb2ycbcr(image)); %转换成Ycbcr颜色空间
image_input = image_input(:,:,1); %取Y通道来进行处理

%定义sobel算子
G_x = [-1,0,1;-2,0,2;-1,0,1];
G_y = [1,2,1;0,0,0;-1,-2,-1];
%定义拉普拉斯算子
L = [-1,-1,-1;-1,8,-1;-1,-1,-1];
image_data = mat2gray(image_input); %图像归一化

for i=2:1:H-1
    for j=2:1:W-1
        image_win = image_data(i-1:i+1, j-1:j+1);
        Gx = sum(sum(image_win .* G_x));
        Gy = sum(sum(image_win .* G_y));
        Sobel_image(i, j) = sqrt(Gx^2 + Gy^2);
        Lapla_image(i, j) = image_data(i,j+1)+image_data(i,j-1)+image_data(i+1,j)+image_data(i-1,j)-4*(image_data(i,j));
    end
end
subplot(131), imshow(image_data), title(\'原图\')
subplot(132), imshow(Sobel_image), title(\'Canny边缘检测\')
subplot(133), imshow(Lapla_image), title(\'拉普拉斯检测\')

代码输出结果如下：

2.4 边缘检测算法总结

Sobel算子检测方法对灰度渐变和噪声较多的图像处理效果较好，sobel算子对边缘定位不是很准确，图像的边缘不止一个像素；当对精度要求不是很高时，是一种较为常用的边缘检测方法。
Canny方法不容易受噪声干扰，能够检测到真正的弱边缘。优点在于，使用两种不同的阈值分别检测强边缘和弱边缘，并且当弱边缘和强边缘相连时，才将弱边缘包含在输出图像中。
Laplacian算子法对噪声比较敏感，所以很少用该算子检测边缘，而是用来判断边缘像素视为与图像的明区还是暗区。拉普拉斯高斯算子是一种二阶导数算子，将在边缘处产生一个陡峭的零交叉, Laplacian算子是各向同性的，能对任何走向的界线和线条进行锐化，无方向性。这是拉普拉斯算子区别于其他算法的最大优点。