边缘检测 Edge detection

前言 

什么是数字图像处理(Digital Image Processing)，先看一下数字图像主要的两个应用领域：

• 改善图示信息以便人们解释；
• 为存储、传输和表示而对图像数据进行处理，以便于机器自动理解

我们可以简单理解为就将一幅原始图像，使用计算机处理为更为我们所能理解或所需要的形式，如图所示，为基于边缘检测的免疫细胞图像自动分割过程示意图

克隆细胞图像自动分割过程示意图

再看一个例子，经典的车牌检测算法，将原始图像进行灰度图转换、边缘检测、形态学腐蚀膨胀等操作，得到车牌区域，随后将车牌区域进行切割(用的是knn，因此识别结果很差)

车牌检测

数字图像处理基础知识与算法

1) 数字图像

数字图像指的是现在的图像都是以二维数字表示，每个像素的灰度值均由一个数字表示，范围为0-255(2^8)

2) 二值图像、灰度图像、彩色图像

二值图像(Binary Image)：图像中每个像素的灰度值仅可取0或1，即不是取黑，就是取白，二值图像可理解为黑白图像

灰度图像(Gray Scale Image)：图像中每个像素可以由0-255的灰度值表示，具体表现为从全黑到全白中间有255个介于中间的灰色值可以取

彩色图像(Color Image)：每幅图像是由三幅灰度图像组合而成，依次表示红绿蓝三通道的灰度值，即我们熟知的RGB，此时彩色图像要视为三维的[height，width, 3]

下面用一张图来感受一下灰度图与彩色图像之间的联系与差别

RGB图像的分解

其中还有一个很重要的公式，即彩色图像转为灰度图的计算公式：

Gray = R×0.299 + G×0.587 + B×0.114 

Gray表示灰度图像，RGB则表示彩色图像的红(red)、绿(green)、蓝(blue)三通道灰度值。

3) 邻接性、连通性

4邻域：假设有一点像素p坐标为(x, y)，则它的4领域是(x + 1, y), (x - 1, y), (x, y + 1), (x, y - 1)

D邻域：假设有一点像素p坐标为(x, y), 则它的D领域是(x + 1, y + 1), ( x + 1, y - 1), (x - 1, y + 1),(x - 1, y - 1)

8邻域：将4领域与D领域的集合取并集，即表示为8邻域

4邻域(左)、 D邻域(中)、 8邻域(右)

4连通：对于在像素点p的4邻域内的像素均与像素点p形成4连通

8连通：对于在像素点p的8邻域内的像素均与像素点p形成8连通

4) 滤波

滤波的目的主要两个：

• 通过滤波来提取图像特征，简化图像所带的信息作为后续其它的图像处理。
• 2为适应图像处理的需求，通过滤波消除图像数字化时所混入的噪声。

其中第一点就是边缘检测中所使用的基本思想，即简化图像信息，使用边缘线代表图像所携带信息

滤波可理解为滤波器(通常为3*3、5*5矩阵)在图像上进行从上到下，从左到右的遍历，计算滤波器与对应像素的值并根据滤波目的进行数值计算返回值到当前像素点，如图所示，蓝色块表示滤波器，对图像进行点积运算并赋值到图像。

具体公式表示为：

 (其中 R5 表示当前像素点， RiGi 表示当前像素与滤波器对应值相乘的值，n为滤波器大小，举例来说如若此滤波器值全为1，则此公式计算的是当前像素点的8连通像素点的平均值，因此滤波完后的图像应表现为模糊的效果，模糊程度取决于滤波器的大小，滤波器大小(size)越大，模糊效果越明显)。

边缘检测(Sobel、Prewitt、Roberts、Canny、Marr-Hildreth)

基本边缘检测算子

边缘检测本质上就是一种滤波算法，区别在于滤波器的选择，滤波的规则是一致的。

为了更好理解边缘检测算子，我们引入梯度(gradient)这一概念，梯度是人工智能(artificial intelligence)非常重要的一个概念，遍布机器学习、深度学习领域，一维函数的一阶微分基本定义为：

图像的滤波一般是基于灰度图进行的，因此图像此时是二维的，看一下二维函数的微分，即偏微分方程：

 由上面的公式我们可以看到，图像梯度即当前所在像素点对于X轴、Y轴的偏导数，所以梯度在图像处理领域可以理解为像素灰度值变化的速度，下面我们举一个简单的例子：

 图中可以看到，100与90之间相差的灰度值为10，即当前像素点在X轴方向上的梯度为10，而其它点均为90，则求导后发现梯度全为0，因此我们可以发现在数字图像处理，因其像素性质的特殊性，微积分在图像处理表现的形式为计算当前像素点沿偏微分方向的差值，所以实际的应用是不需要用到求导的，只需进行简单的加减运算。

而另一个概念梯度的模则表示f(x, y)，在其最大变化率方向上的单位距离所增加的量，即：

在了解完梯度的概念之后呢，下面我们先介绍一下几种基本边缘检测滤波器: Sobel、Prewitt、Roberts算子。

左边：Roberts算子；右边：Prewitt算子。

Sobel算子。

滤波算子作用过程

以Sobel为例，其中 Sx,Sy 分别表示对于X轴、Y轴的边缘检测算子，从 Sx 算子结构可以很清楚发现，这个滤波器是计算当前像素点右边与左边8连通像素灰度值的差值，我们先通过一维的概念来理解一下：

如现在有一个一维数组长度为10，值为：

[ 8, 6, 2, 4, 9, 1, 3, 5, 10, 6 ]

此时我们的一维边缘检测算子则表现为[ -1, 0, 1 ]，现在我们把边缘检测算子放在数组上面进行点积(即对应点相乘之后的和)，得到结果为：

[ 6, -6, -2, 7, -3, -6, 4, 7, 1, -10]

可以发现得到的值出现了负数，但我们之前的定义则声明了像素灰度值定义域为0-255范围内，因此这里一般的操作为将负数截断到0-255以内或者直接取绝对值，因此现在我们得到的是

[ 6, 6, 2, 7, 3, 6, 4, 7, 1, 10]

其中数字的大小则表示了当前像素点梯度的模大小，即灰度变化的速度有多大，值越大，我们一定程度上就可以确信当前点为我们所要找的边缘点，通过一维的例子我们可以更好理解二维的边缘检测思想，即沿着X轴、Y轴进行两次滤波操作，得到的结果进行平方求和加根号的操作得出当前像素点的图像梯度，我们来通过一张图理解一下这个过程：

原图像、沿X轴梯度图像、沿Y轴梯度图像、梯度图像的可视化.

图中(a)为原始的灰度图像，(b)和(c)为使用图3-3中Sobel算子的 、Sx、Sy 两种形式分别对原始图像进行的滤波结果，即表示为分别沿X、Y轴的梯度图，最后将两个图融合在一起则得到了我们所需的梯度图像(d)，在给大家一张图来帮助理解Sobel算法。

现在我们已经大致了解了边缘检测的基本思想了，直接用基本的边缘算子如Sobel求得的边缘图存在很多问题，如噪声污染没有被排除、边缘线太过于粗宽等，因此我们接下来要介绍两个先进的边缘检测算子：Canny算子和Marr-Hildreth算子。

较为先进的边缘检测算子

1）Canny算子

Canny算子是澳洲计算机科学家约翰·坎尼（John F. Canny）于1986年开发出来的一个多级

边缘检测算法，其目标是找到一个最优的边缘，其最优边缘的定义是：

• 好的检测 --算法能够尽可能多地标示出图像中的实际边缘。
• 好的定位 --标识出的边缘要与实际图像中的实际边缘尽可能接近。
• 最小响应 --图像中的边缘只能标识一次，并且可能存在的图像噪声不应该标识为边缘。

所以接下来我们来介绍一下目前流行的Canny算法的具体步骤：

（1）高斯(Gaussian)滤波

高斯滤波目前是最为流行的去噪滤波算法，高斯与我们学的概率论中正态分布中正态一词指的是同一个意思，其原理为根据待滤波的像素点及其邻域点的灰度值按照高斯公式生成的参数规则进行加权平均，这样可以有效滤去理想图像中叠加的高频噪声(noise)

二维高斯公式为：

 常见的高斯滤波器有：

 其实高斯滤波器很像一个金字塔结构，其滤波器的值大小我们可以理解为权重(weight)，值越大对应的像素点权重越大，分量也就越大，因此从高斯滤波器我们可以看出对应当前像素点，距离越远权重越小，对灰度值的贡献也就越小。

让我们举个例子来理解一下高斯滤波，如上图左边的高斯滤波器，其中心点4我们可以把它看成是'主人公'，其周围的点看成是'邻居'，噪声我们把它看成是'坏人'，现在我们假设这9个人里面，有一个人是'坏人'，我们也知道坏人是肯定会说自己是好人的，但要是我们有投票机制决定一个人是否为'坏人'呢，其中权重(weight)则对应每个人说话的分量，投票机制就为我们所说的加权平均策略，现在我们可以很直观地发现，其实高斯滤波就是一个会考虑其周围像素点的滤波器，即使当前点位为噪声点，高斯滤波器也会通过周围点的灰度值来制约噪声的影响，生成高斯滤波器与滤波的代码如下：

sigma=1;               %高斯标准差

% 根据高斯标准差计算滤波器长度
filterExtent = ceil(4*sigma);
x = -filterExtent:filterExtent;

% 生成一维高斯核
c = 1/(sqrt(2*pi)*sigma);
gaussKernel = c * exp(-(x.^2)/(2*sigma^2));

% 标准化
gaussKernel = gaussKernel/sum(gaussKernel);

% 对图像进行高斯滤波平滑
aSmooth=imfilter(a,gaussKernel,'conv','replicate');   % 沿着X轴卷积
aSmooth=imfilter(aSmooth,gaussKernel','conv','replicate');  % 沿着Y轴卷积

（其中gaussKernel'表示对gaussKernel进行转置）

（2）计算梯度图像与角度图像

计算梯度图像我们刚才基本也有理解了一下，无非就是用各种边缘检测算子进行梯度的检测，但Canny中使用的梯度检测算子有点高级，是使用高斯滤波器进行梯度计算得到的滤波器，得到的结果也类似于Sobel算子，即距离中心点越近的像素点权重越大，代码如下：

% 数值梯度函数(Gaussian核的生成1-D导数)
derivGaussKernel = gradient(gaussKernel);

% 标准化
negVals = derivGaussKernel < 0;
posVals = derivGaussKernel > 0;
derivGaussKernel(posVals) = derivGaussKernel(posVals)/sum(derivGaussKernel(posVals));
derivGaussKernel(negVals) = derivGaussKernel(negVals)/abs(sum(derivGaussKernel(negVals)));

% 计算梯度
dx = imfilter(aSmooth, derivGaussKernel, 'conv','replicate');
dy = imfilter(aSmooth, derivGaussKernel', 'conv','replicate');


mag = hypot(dx,dy); 
magmax = max(mag(:));
if magmax>0
    magGrad = mag / magmax;   % 梯度标准化
end

角度图像的计算则较为简单，其作用为非极大值抑制的方向提供指导，公式如下：

（3）对梯度图像进行非极大值抑制

从上一步得到的梯度图像存在边缘粗宽、弱边缘干扰等众多问题，现在我们可以使用非极大值抑制来寻找像素点局部最大值，将非极大值所对应的灰度值置0，这样可以剔除一大部分非边缘的像素点。

如下图所示，C表示为当前非极大值抑制的点，g1-4为它的8连通邻域点，图中蓝色线段表示上一步计算得到的角度图像C点的值，即梯度方向，第一步先判断C灰度值在8值邻域内是否最大，如是则继续检查图中梯度方向交点dTmp1,dTmp2值是否大于C，如C点大于dTmp1,dTmp2点的灰度值，则认定C点为极大值点，置为1，因此最后生成的图像应为一副二值图像，边缘理想状态下都为单像素边缘。

非极大值抑制

 （其中需要注意的是梯度方向交点并不一定落在8领域所在8个点的位置，因此dTmp1和dTmp2实际应用中是使用相邻两个点的双线性插值所形成的灰度值）。

最后在上一张图帮助大家理解，如下图所示，其中梯度方向均为垂直向上，经过非极大值抑制后取梯度方向上最大值为边缘点，形成细且准确的单像素边缘。

 （4）使用双阈值进行边缘连接

经过以上三步之后得到的边缘质量已经很高了，但还是存在很多伪边缘，因此Canny算法中所采用的算法为双阈值法，具体思路为选取两个阈值，将小于低阈值的点认为是假边缘置0，将大于高阈值的点认为是强边缘置1，介于中间的像素点需进行进一步的检查。

根据高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，该算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像闭合，具体代码为：

function nedge=connect1(nedge,y,x,low,high,magGrad)       %种子定位后的连通分析
    neighbour=[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];  %八连通搜寻
    [m n]=size(nedge);

    for k=1:8
        yy=fix(y+neighbour(k,1));
        xx=fix(x+neighbour(k,2));


        if yy>=1 &&yy<=m &&xx>=1 && xx<=n  
            if magGrad(yy,xx)>=low & nedge(yy,xx)~=255 & magGrad(yy,xx)<high  
                nedge(yy,xx)=255;
                %disp('check check');
                %nedge=connect1(nedge,yy,xx,low,high,magGrad);
            end
        end
    end
end

但由于寻找弱边缘点的计算代价过大，因为使用的是递归思维，且所找寻到的弱边缘点为数不多，因此实际应用中常常舍去这一步骤，取而代之的是基于形态学的边缘细化操作，具体思想我们以后还会提到，具体代码为：

H = bwmorph(H, 'thin', 1); 

至此，我们已经深度理解了Canny算法的思想与实现手段，实际应用中Canny一般是边缘检测的首选项，其算法思想也非常值得我们学习，接下来我们在简单介绍基于二阶导数法的Marr-Hildreth边缘检测算子。

2）Marr-Hildreth算子

在学习Marr-Hildreth算子之前我们先来理解一下为什么要用二阶导数法

如图 3-8所示，左边表示的是一副灰度图像，从左到右从黑色(0)慢慢变为白色(255)，现在我们来看它的水平灰度剖面，灰度值从小到大平稳上升，其一阶导数表示为在上升区域为不变的值，其中得到的信息是图像灰度值是平稳过渡的，即梯度值相等，接下来在将其求二次导数，得到的图像为在开始过渡的起点为正数，其值为一阶导数在此点的梯度值，结束点也和起点一样，现在重点来了，我们将这两点连起来得到一个与X轴的交叉点，这一点就是我们所认为的边缘点，这就是二阶导数应用在边缘检测领域的奇妙之处（第一次接触的时候觉得巨神奇）

看一下marr和hildreth的证明结论：

• 灰度变化与图像尺寸无关，因此他们的检测要求使用不同尺寸的算子。
• 灰度的突然变换会在一阶导数中引起波峰或波谷，在二阶导数中等效引起零交叉。

 

 一条边缘的两个附加性质：

• 对图像中的每条边缘，二阶导数生成两个值
•   二阶导数的零交叉点可用于定位粗边缘的中心

总结：二阶导数提取边缘往往产生双边缘。而 通过求零交叉点可以确定边缘的中心，从而避免了产生双边缘的不方便。

来看看思路。

(1) 高斯滤波

基本所有边缘检测算法前面都会加一个高斯滤波来去除高频噪声。

(2) 计算拉普拉斯(Laplacian)二阶导

Marr-Hildreth证明，最满足图像处理需求的算子是滤波器拉普拉斯高斯(Log)算子，具体原理不多阐述，这里来看一下它的公式：

 由其生成的拉普拉斯滤波器也被称为墨西哥草帽算子，因为其中间一般为较大的正数，8邻域连通点为较小的负数值，常用的滤波器如图所示：

拉普拉斯滤波器

 之后就是使用拉普拉斯滤波器进行图像的滤波操作，得到待计算图像。

(3) 计算零交叉(Zero crossing)

零交叉的实现较为简单，由于零交叉点意味着至少两个相邻的像素点的像素值异号，一共有四种需要检测的情况：左右，上下，两个对角，其中如果滤波后的图像g(x, y)的任意像素p处的四种情况其中一组的差值的绝对值超过了设定的阈值，则我们可以称p为一个零交叉像素，示例如下：

 此为Marr-Hildreth其中一小部分，检测[ - +]这一情况是否满足，其中thresh为提到的阈值。

到这里学习完两种最为流行且经典的先进边缘检测算法与思想，接下来说的是一些经验与算法的选择参考。

补充

1. 滤波器的大小应该是奇数，这样才有一个中心点可进行赋值操作，常见的滤波器或卷积核(Conv kernel)有3*3、5*5等，因此也有了半径的概念，例如5*5的卷积核的半径为2。
2. 滤波器中所有元素之和应为0，这一限制条件是保证滤波前后图像总体灰度值不变。
3. Roberts算子、Sobel算子、Prewitt算子运算速率高，对噪声也有一定抑制作用，但检测出的边缘质量不高，如边缘较粗、定位不准、间断点多。
4. Canny算子不容易受噪声干扰，得到的边缘精细且准确，缺点就是运算代价较高，运行于实时图像处理较困难，适用于高精度要求的应用。
5. Marr-Hildreth算子边缘检测效果相对较优，但对于噪声比较敏感（因其二阶运算的性质）。

Canny边缘检测算法完整的Matlab代码实现

I=imread('cameraman.tif');
%I=rgb2gray(I);
figure(1);subplot(121);imshow(I);xlabel('原图像');
[m n]=size(I);
a=double(I);
sigma=1;               %高斯标准差
%highThresh=0.0625;     %上阈值
%lowThresh=0.0250 ;    %下阈值
  
  
%=======================高斯滤波&梯度计算=======================

%%%%%%%%%%%%%%%%%%%%%%%%%Old%%%%%%%%%%%%%%%%%%%%%%%%%%%
%pw = 1:30; 
%ssq = sigma^2;
%width = find(exp(-(pw.*pw)/(2*ssq))>0.0001,1,'last');
%if isempty(width)
%    width = 1;  
%end
    
%t = (-width:width);
%gauss = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);     % 一维高斯滤波器
    

%[x,y]=meshgrid(-width:width,-width:width);
%gauss2=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);   %二维高斯滤波器
    
%对图像进行高斯滤波平滑
%aSmooth=imfilter(a,gauss,'conv','replicate');   % 沿着X轴卷积
%aSmooth=imfilter(aSmooth,gauss','conv','replicate'); % 沿着Y轴卷积
    
%使用二维高斯对图像进行卷积
%dx = imfilter(aSmooth, gauss2, 'conv','replicate');
%dy = imfilter(aSmooth, gauss2', 'conv','replicate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% 根据高斯标准差计算滤波器长度
filterExtent = ceil(4*sigma);
x = -filterExtent:filterExtent;

% 生成一维高斯核
c = 1/(sqrt(2*pi)*sigma);
gaussKernel = c * exp(-(x.^2)/(2*sigma^2));

% 标准化
gaussKernel = gaussKernel/sum(gaussKernel);

% 数值梯度函数(Gaussian核的生成1-D导数)
derivGaussKernel = gradient(gaussKernel);

% 标准化
negVals = derivGaussKernel < 0;
posVals = derivGaussKernel > 0;
derivGaussKernel(posVals) = derivGaussKernel(posVals)/sum(derivGaussKernel(posVals));
derivGaussKernel(negVals) = derivGaussKernel(negVals)/abs(sum(derivGaussKernel(negVals)));

    
% 对图像进行高斯滤波平滑
aSmooth=imfilter(a,gaussKernel,'conv','replicate');   % 沿着X轴卷积
aSmooth=imfilter(aSmooth,gaussKernel','conv','replicate');  % 沿着Y轴卷积
%hv=fspecial('sobel');

% 计算梯度
dx = imfilter(aSmooth, derivGaussKernel, 'conv','replicate');
dy = imfilter(aSmooth, derivGaussKernel', 'conv','replicate');


mag = hypot(dx,dy); 
magmax = max(mag(:));
if magmax>0
    magGrad = mag / magmax;   % 梯度标准化
end

% 阈值选择
%PercentOfPixelsNotEdges = 0.7; 
counts=imhist(magGrad, 64);
highThresh = find(cumsum(counts) > 0.7*m*n, 1 ,'first' ) / 64;
lowThresh = 0.4*highThresh;

%figure(8);imshow(magGrad);
%%========================高斯滤波========================================
%w=fspecial('gaussian',[5 5]);
%img=imfilter(img,w,'replicate');
%figure;
%imshow(uint8(img))

%%===================sobel边缘检测=======================================
%hv=fspecial('sobel');
%dx=imfilter(img,hv,'replicate');      %求横边缘
%hh=hv';
%dy=imfilter(img,hh,'replicate');      %求竖边缘
%img=sqrt(dx.^2+dy.^2); 


%    magmax = max(img(:));%      （阈值选择归一化）
%    if magmax > 0
%        magGrad = img / magmax;
%    end
%figure;
%imshow(uint8(img));

I = thinAndThreshold(dx, dy, magGrad, lowThresh, highThresh);


%disp(lowThresh);
subplot(122);imshow(I);xlabel('canny边缘检测');
disp("高阈值TL: "+highThresh);
disp("低阈值TH: "+lowThresh);


%========================非极大值抑制和边缘连接=======================================
function H = thinAndThreshold(dx, dy, magGrad, lowThresh, highThresh)

E = cannyFindLocalMaxima(dx,dy,magGrad,lowThresh);  %非极大值抑制

if ~isempty(E)
    [rstrong,cstrong] = find(magGrad>highThresh & E);

    if ~isempty(rstrong) 
        H = bwselect(E, cstrong, rstrong, 8);   % 选定强边缘8连通目标
       % figure(2);imshow(H);
        
       % set(0,'RecursionLimit',1000);      %弱边缘连通（无太大作用，且运算时间过长）
       % [xstrong ystrong]=find(magGrad>highThresh & E);
       % for i=1:numel(xstrong);
       %     H = connect1(H,xstrong(i),ystrong(i),lowThresh,highThresh,magGrad);
       % end
       %  figure(3);imshow(H);
        
        H = bwmorph(H, 'thin', 1);      % 边缘细化
    else
        H = false(size(E));
    end
else
    H = false(size(E));
end
end

%========================弱边缘连接=======================================

function nedge=connect1(nedge,y,x,low,high,magGrad)       %种子定位后的连通分析
    neighbour=[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];  %八连通搜寻
    [m n]=size(nedge);

    for k=1:8
        yy=fix(y+neighbour(k,1));
        xx=fix(x+neighbour(k,2));


        if yy>=1 &&yy<=m &&xx>=1 && xx<=n  
            if magGrad(yy,xx)>=low & nedge(yy,xx)~=255 & magGrad(yy,xx)<high  
                nedge(yy,xx)=255;
                %disp('check check');
                %nedge=connect1(nedge,yy,xx,low,high,magGrad);
            end
        end
    end
end

左图：原图；右图：边缘图。

参考博客

数字图像处理:边缘检测(Edge detection) - 知乎 (zhihu.com)