文章目录

Roberts、Prewitt、Sobel、Laplacian、LoG 和 Canny 边缘检测算子(MATLAB自写函数实现)
- 1理论
- - 1.1 知识引入
  - - 1.1.1 图像边缘边缘[1]
    - 1.1.2 一阶导数[2]
    - 1.1.3 二阶导数[2]
    - 1.1.4 一阶导数与二阶导数区别
  - 1.2 边缘检测算子
  - - 1.2.1 Roberts[1][2]
    - 1.2.2Prewitt[1][2]
    - 1.2.3 Sobel [1][2]
    - 1.2.4Laplacian[1]
    - 1.2.5 LoG [4]
    - 1.2.6 Canny [5]
- 2 上代码！
- - 2.1 边界处理函数
  - 2.2 Canny算法函数
  - 2.3 边缘检测函数
- 3 demo结果

Roberts、Prewitt、Sobel、Laplacian、LoG 和 Canny 边缘检测算子(MATLAB自写函数实现)

1理论

1.1 知识引入

1.1.1 图像边缘边缘[1]

图像边缘是图像最基本的特征，所谓边缘(Edge) 是指图像局部特性的不连续性。灰度或结构等信息的突变处称之为边缘。例如，灰度级的突变、颜色的突变,、纹理结构的突变等。边缘是一个区域的结束，也是另一个区域的开始，利用该特征可以提取图像边缘。

图(a)是一个理想的边缘所具备的特性。每个灰度级跃变到一个垂直的台阶上。而实际上，在图像采集系统的性能、采样率和获取图像的照明条件等因素的影响，得到的边缘往往是模糊的，边缘被模拟成具有“斜坡面”的剖面，如图(b)所示，在这个模型中，模糊的边缘变得“宽”了，而清晰的边缘变得“窄”了。

图像的边缘有方向和幅度两种属性。边缘通常可以通过一阶导数或二阶导数检测得到。一阶导数是以最大值作为对应的边缘的位置，而二阶导数则以过零点作为对应边缘的位置。

1.1.2 一阶导数[2]

在一维连续数集上有函数 f ( x ) f(x) f(x),我们可以通过求导获得该函数在任一点的斜率，根据导数的定义有(以下 Δ x \Delta x Δx， Δ y \Delta y Δy都趋于零)：
f ′ ( x ) = f ( x + Δ x ) − f ( x ) Δ x f'(x)=\frac{f(x+\Delta x)-f(x)}{\Delta x} f′(x)=Δxf(x+Δx)−f(x)
在二维连续数集上有函数 f ( x , y ) f(x,y) f(x,y),我们也可以通过求导获得该函数在 x x x和 y y y分量的偏导数，根据定义有：
g x = ∂ f ( x , y ) ∂ x = f ( x + Δ x , y ) − f ( x , y ) Δ x g_x=\frac{\partial f(x,y)}{\partial x}=\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x} gx=∂x∂f(x,y)=Δxf(x+Δx,y)−f(x,y)
g y = ∂ f ( x , y ) ∂ y = f ( x , y + Δ y ) − f ( x , y ) Δ y g_y=\frac{\partial f(x,y)}{\partial y}=\frac{f(x,y+\Delta y)-f(x,y)}{\Delta y} gy=∂y∂f(x,y)=Δyf(x,y+Δy)−f(x,y)

对于图像来说，是一个二维的离散型数集，通过推广二维连续型求函数偏导的方法，来求得图像的偏导数：

g x = ∂ f ( x , y ) ∂ x = f ( x + 1 , y ) − f ( x , y ) g_x=\frac{\partial f(x,y)}{\partial x}=f(x+1,y)-f(x,y) gx=∂x∂f(x,y)=f(x+1,y)−f(x,y)
g y = ∂ f ( x , y ) ∂ x = f ( x , y + 1 ) − f ( x , y ) g_y=\frac{\partial f(x,y)}{\partial x}=f(x,y+1)-f(x,y) gy=∂x∂f(x,y)=f(x,y+1)−f(x,y)
若以图像左上角点为图像坐标原点，原点向右为 x x x轴，原点向下为 y y y轴，则
g x g_x gx对应的模板为：

g y g_y gy对应的模板为：

梯度对应一阶导数﹐梯度算子是一阶导数算子。对一个连续函数 f ( x , y ) f(x,y) f(x,y)，它在位置 ( x , y ) (x,y) (x,y)的梯度可表示为一个向量(两个分量分别是沿 X X X和 Y Y Y方向的一阶导数)
∇ f ( x , y ) = [ g x , g y ] T = [ ∂ f ∂ x ， ∂ f ∂ y ] \nabla f(x,y)=[g_x,g_y]^T=[\frac{\partial f}{\partial x}，\frac{\partial f}{\partial y}] ∇f(x,y)=[gx,gy]T=[∂x∂f，∂y∂f]

这个向量的幅值(也常直接简称为梯度)和方向角分别为
幅值： M ( x , y ） = ∣ ∣ ∇ f ∣ ∣ 2 = [ g x 2 + g y 2 ] 1 / 2 幅值：M(x,y）=||\nabla f||_2=[g_x^2+g_y^2]^{1/2} 幅值：M(x,y）=∣∣∇f∣∣2=[gx2+gy2]1/2
方向角： ϕ ( x , y ) = a r c t a n ( g y / g x ) 方向角：\phi(x,y)=arctan(g_y/g_x) 方向角：ϕ(x,y)=arctan(gy/gx)
由于幅值计算的是2范数涉及平方和开方运算,计算量会比较大。在实用中为了计算简便,常采用其他方法组合两个模板的输出。有一种简单的方法是求1范数的,即
∣ ∣ ∇ f ∣ ∣ 1 = ∣ g x ∣ + ∣ g y ∣ ||\nabla f||_1=|g_x|+|g_y| ∣∣∇f∣∣1=∣gx∣+∣gy∣

另一种简单的方法是求 ∞ \infty ∞范数,即
∣ ∣ ∇ f ∣ ∣ ∞ = m a x { ∣ g x ∣ , ∣ g y ∣ } ||\nabla f||_{\infty}=max\{|g_x|,|g_y|\} ∣∣∇f∣∣∞=max{∣gx∣,∣gy∣}

上面是图像的垂直和水平梯度，但我们有时候也需要对角线方向的梯度，定义如下：
g x = ∂ f ( x , y ) ∂ x = f ( x + 1 , y + 1 ) − f ( x , y ) g_x=\frac{\partial f(x,y)}{\partial x}=f(x+1,y+1)-f(x,y) gx=∂x∂f(x,y)=f(x+1,y+1)−f(x,y)
g y = ∂ f ( x , y ) ∂ x = f ( x + 1 , y ) − f ( x , y + 1 ) g_y=\frac{\partial f(x,y)}{\partial x}=f(x+1,y)-f(x,y+1) gy=∂x∂f(x,y)=f(x+1,y)−f(x,y+1)

1.1.3 二阶导数[2]

一阶方向导数：
∂ f ( x , y ) ∂ x = f ( x + 1 , y ) − f ( x , y ) \frac{\partial f(x,y)}{\partial x}=f(x+1,y)-f(x,y) ∂x∂f(x,y)=f(x+1,y)−f(x,y)
∂ f ( x , y ) ∂ y = f ( x , y + 1 ) − f ( x , y ) \frac{\partial f(x,y)}{\partial y}=f(x,y+1)-f(x,y) ∂y∂f(x,y)=f(x,y+1)−f(x,y)
二阶方向偏导：

∂ 2 f ( x , y ) ∂ x 2 = f ′ ( x + 1 , y ) − f ′ ( x , y ) = f ( x + 2 , y ) − f ( x + 1 , y ) − ( f ( x + 1 , y ) − f ( x , y ) ) = f ( x + 2 , y ) − 2 f ( x + 1 , y ) + f ( x , y ) \frac{\partial^2f(x,y)}{\partial x^2}=f'(x+1,y)-f'(x,y)\\=f(x+2,y)-f(x+1,y)-(f(x+1,y)-f(x,y))\\=f(x+2,y)-2f(x+1,y)+f(x,y) ∂x2∂2f(x,y)=f′(x+1,y)−f′(x,y)=f(x+2,y)−f(x+1,y)−(f(x+1,y)−f(x,y))=f(x+2,y)−2f(x+1,y)+f(x,y)

∂ 2 f ( x , y ) ∂ y 2 = f ′ ( x , y + 1 ) − f ′ ( x , y ) = f ( x , y + 2 ) − f ( x , y + 1 ) − ( f ( x , y + 1 ) − f ( x , y ) ) = f ( x , y + 2 ) − 2 f ( x , y + 1 ) + f ( x , y ) \frac{\partial^2f(x,y)}{\partial y^2}=f'(x,y+1)-f'(x,y)\\=f(x,y+2)-f(x,y+1)-(f(x,y+1)-f(x,y))\\=f(x,y+2)-2f(x,y+1)+f(x,y) ∂y2∂2f(x,y)=f′(x,y+1)−f′(x,y)=f(x,y+2)−f(x,y+1)−(f(x,y+1)−f(x,y))=f(x,y+2)−2f(x,y+1)+f(x,y)

将上式中的变量减1后，得到：
∂ 2 f ( x , y ) ∂ x 2 = f ( x + 1 , y ) − 2 f ( x , y ) + f ( x − 1 , y ) \frac{\partial^2f(x,y)}{\partial x^2}=f(x+1,y)-2f(x,y)+f(x-1,y) ∂x2∂2f(x,y)=f(x+1,y)−2f(x,y)+f(x−1,y)

∂ 2 f ( x , y ) ∂ x 2 = f ( x , y + 1 ) − 2 f ( x , y ) + f ( x , y − 1 ) \frac{\partial^2f(x,y)}{\partial x^2}=f(x,y+1)-2f(x,y)+f(x,y-1) ∂x2∂2f(x,y)=f(x,y+1)−2f(x,y)+f(x,y−1)

二阶导数，其定义如下：

∇ 2 f ( x , y ) = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 = [ f ( x + 1 , y ) − 2 f ( x ) + f ( x − 1 ) ] + [ f ( x , y + 1 ) − 2 f ( x ) + f ( x , y − 1 ) ] = f ( x + 1 , y ) + f ( x − 1 , y ) + f ( x , y + 1 ) + f ( x , y − 1 ) − 4 f ( x , y ) \nabla^2f(x,y)=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}\\=[f(x+1,y)-2f(x)+f(x-1)]+[f(x,y+1)-2f(x)+f(x,y-1)]\\=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y) ∇2f(x,y)=∂x2∂2f+∂y2∂2f=[f(x+1,y)−2f(x)+f(x−1)]+[f(x,y+1)−2f(x)+f(x,y−1)]=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)
其对应的模板：
( 1 0 1 0 − 4 0 1 0 1 ) \left(\begin{array}{ccc}1 & 0 & 1 \\0 & -4 & 0 \\1 & 0 & 1\end{array}\right) ⎝⎛1010−40101⎠⎞

1.1.4 一阶导数与二阶导数区别

设上图最上面部分的灰度剖面图对应于图像中的一条具有代表性的水平像素线，其中包括了

灰度较缓变化的斜坡（软边缘）
孤立点（很可能为噪声）
细线（细节）
以及灰度跳变的阶梯（硬边缘）

为了简单起见，考虑图像中只有8个灰度级（0到7）的情况。图中下面的一行给出了这条像素线中各个像素的灰度值，由此计算出的一阶微分和二阶微分在图中的第3行和第4行中给出。由于这里的像素线在图像中是水平分布的，因此微分方程也是一维的，一维情况下的一阶微分：
∂ f ∂ x = f ( x + 1 ) − f ( x ) \frac{\partial f}{\partial x}=f(x+1)-f(x) ∂x∂f=f(x+1)−f(x)

和一维情况下的二阶微分：
∂ 2 f ∂ x 2 = f ( x + 1 ) + f ( x − 1 ) − 2 f ( x ) \frac{\partial^{2} f}{\partial x^{2}}=f(x+1)+f(x-1)-2 f(x) ∂x2∂2f=f(x+1)+f(x−1)−2f(x)

通过分析这个典型的灰度变化模型，就可以很好地比较噪声点、细节以及边缘的一阶和二阶微分结果。

首先注意到沿着整个斜坡（软边缘），一阶微分都具有非0响应，并且当这种斜坡的灰度过渡近似线时，对应于变化率的一阶微分的响应为恒定值（这里为-1）；而二阶微分的非0响应则只出现在斜坡的起始和终点处，在灰度变化率恒定的斜面上二阶微分值为0，这就是拉普拉斯锐化图像周围出现双边缘的原因。由此得出结论，对于图像中的软边缘，一阶微分通常产生较粗的边缘，而二阶微分则细得多。

再来看孤立噪声点，注意到二阶微分对于噪声点的响应较一阶微分要强很多，这也就是拉普拉斯锐化图像中出现一些零星的高响应的原因，当然二阶微分的这一性质是我们所不希望的。

细线常常对应于图像中的细节，二阶微分对细线的较强响应说明了二阶微分对于细节增强的优越性。

最后，一、二阶微分对于灰度阶梯有着相同的响应，只是在二阶微分中有一个从正到负的过渡，这一性质用于边缘检测。

前面介绍了一些图像边缘检测相关知识，接下来介绍常用的几个边缘检测算子。

1.2 边缘检测算子

1.2.1 Roberts[1][2]

Roberts算子又称为交叉微分算法，它是基于交叉差分的梯度算法，通过局部差分计算检测边缘线条。Roberts算子的模板分为水平方向和垂直方向，能较好的增强正负45度的图像边缘。

下面给出Roberts算子的模板，在像素点P5处 x x x 和 y y y 方向上的梯度大小 g x g_{x} gx和 g y g_{y} gy 分别计算为：

g x = ∂ f ∂ x = P 9 − P 5 g_{x}=\frac{\partial f}{\partial x}=\mathrm{P} 9-\mathrm{P} 5 gx=∂x∂f=P9−P5
g y = ∂ f ∂ y = P 8 − P 6 g_{y}=\frac{\partial f}{\partial y}=\mathrm{P} 8-\mathrm{P} 6 gy=∂y∂f=P8−P6

1.2.2Prewitt[1][2]

Prewitt算子是一种图像边缘检测的微分算子，其原理是利用特定区域内像素灰度值产生的差分实现边缘检测。由于Prewitt算子采用 3 × 3 3\times3 3×3 模板对区域内的像素值进行计算，其边缘检测结果在水平方向和垂直方向均比Robert算子更加明显。
下面给出Prewitt算子的模板，在像素点P5处 x x x 和 y y y 方向上的梯度大小 g x g_{x} gx 和 g y g_{y} gy分别计算为：

g x = ∂ f ∂ x = ( P 7 + P 8 + P 9 ) − ( P 1 + P 2 + P 3 ) g_{x}=\frac{\partial f}{\partial x}=(\mathrm{P} 7+\mathrm{P} 8+\mathrm{P} 9)-(\mathrm{P} 1+\mathrm{P} 2+\mathrm{P} 3) gx=∂x∂f=(P7+P8+P9)−(P1+P2+P3)

g y = ∂ f ∂ y = ( P 3 + P 6 + P 9 ) − ( P 1 + P 4 + P 7 ) g_{y}=\frac{\partial f}{\partial y}=(\mathrm{P} 3+\mathrm{P} 6+\mathrm{P} 9)-(\mathrm{P} 1+\mathrm{P} 4+\mathrm{P} 7) gy=∂y∂f=(P3+P6+P9)−(P1+P4+P7)

也可以沿对角计算：

1.2.3 Sobel [1][2]

Sobel算子在Prewitt算子的基础上增加了权重的概念，认为相邻点的距离远近对当前像素点的影响是不同的，距离越近的像素点对应当前像素的影响越大，从而实现图像锐化并突出边缘轮廓。

g x = ∂ f ∂ x = ( P 7 + 2 P 8 + P 9 ) − ( P 1 + 2 P 2 + P 3 ) g_{x}=\frac{\partial f}{\partial x}=(\mathrm{P} 7+2 \mathrm{P} 8+\mathrm{P} 9)-(\mathrm{P} 1+2 \mathrm{P} 2+\mathrm{P} 3) gx=∂x∂f=(P7+2P8+P9)−(P1+2P2+P3)

g y = ∂ f ∂ y = ( P 3 + 2 P 6 + P 9 ) − ( P 1 + 2 P 4 + P 7 ) g_{y}=\frac{\partial f}{\partial y}=(\mathrm{P} 3+2 \mathrm{P} 6+\mathrm{P} 9)-(\mathrm{P} 1+2 \mathrm{P} 4+\mathrm{P} 7) gy=∂y∂f=(P3+2P6+P9)−(P1+2P4+P7)

也可以沿对角计算：

1.2.4Laplacian[1]

拉普拉斯(Laplacian) 算子是 n 维欧几里德空间中的一个二阶微分算子，分为四邻域和八邻域，四邻域是对邻域中心像素的四个方向求梯度，八邻域是对八个方向求梯度。
其中，Laplacian算子四邻域模板如下所示：
H = [ 0 − 1 0 − 1 4 − 1 0 − 1 0 ] \mathrm{H}=\left[\begin{array}{ccc}0 & -1 & 0 \\-1 & 4 & -1 \\0 & -1 & 0\end{array}\right] H=⎣⎡0−10−14−10−10⎦⎤

Laplacian算子八邻域模板如下所示：
H = [ − 1 − 1 − 1 − 1 8 − 1 − 1 − 1 − 1 ] \mathrm{H}=\left[\begin{array}{ccc}-1 & -1 & -1 \\-1 & 8 & -1 \\-1 & -1 & -1\end{array}\right] H=⎣⎡−1−1−1−18−1−1−1−1⎦⎤
注意：模板系数之和为0，是为了当遇到没有边缘的区域时，模板响应为0。

1.2.5 LoG [4]

由于拉普拉斯算子是一个二阶导数，对噪声具有无法接受的敏感性，而且其幅值会产生双边缘，另外，边缘方向的不可检测性也是拉普拉斯算子的缺点之一，因此，一般不以其原始形式用于边缘检测。

为了弥补拉普拉斯算子与生俱来的缺陷，美国学者Marr提出了一种算法，在运用拉普拉斯算子之前一般先进行高斯低通滤波，可表示为：
∇ 2 [ G ( x , y ) ∗ f ( x , y ) ] \nabla^{2}[G(x, y) * f(x, y)] ∇2[G(x,y)∗f(x,y)]

其中 f ( x , y ) f(x, y) f(x,y)为图像， G ( x , y ) G(x, y) G(x,y)为高斯函数，表示为：
G ( x , y ) = 1 2 π σ 2 exp ⁡ ( − x 2 + y 2 2 σ 2 ) G(x, y)=\frac{1}{2 \pi \sigma^{2}} \exp \left(-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right) G(x,y)=2πσ21exp(−2σ2x2+y2)

其中， σ σ σ是标准差。用高斯函数卷积模糊一幅图像，图像模糊的程度是由 σ σ σ决定的。

由于在线性系统中卷积与微分的次序可以交换，可得下式：
∇ 2 [ G ( x , y ) ∗ f ( x , y ) ] = ∇ 2 G ( x , y ) ∗ f ( x , y ) \nabla^{2}[G(x, y) * f(x, y)]=\nabla^{2} G(x, y) * f(x, y) ∇2[G(x,y)∗f(x,y)]=∇2G(x,y)∗f(x,y)

上式说明了可以先对高斯算子进行微分运算，然后再与图像 f ( x , y ) f(x,y) f(x,y)卷积，其效果等价于在运用拉普拉斯算子之前首先进行高斯低通滤波。

计算式 G ( x , y ) G(x,y) G(x,y)的二阶偏导，如下式所示。
∂ 2 G ( x , y ) ∂ x 2 = 1 2 π σ 4 [ x 2 σ 2 − 1 ] exp ⁡ ( − x 2 + y 2 2 σ 2 ) \frac{\partial^{2} G(x, y)}{\partial x^{2}}=\frac{1}{2 \pi \sigma^{4}}\left[\frac{x^{2}}{\sigma^{2}}-1\right] \exp \left(-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right) ∂x2∂2G(x,y)=2πσ41[σ2x2−1]exp(−2σ2x2+y2)

∂ 2 G ( x , y ) ∂ y 2 = 1 2 π σ 4 [ y 2 σ 2 − 1 ] exp ⁡ ( − x 2 + y 2 2 σ 2 ) \frac{\partial^{2} G(x, y)}{\partial y^{2}}=\frac{1}{2 \pi \sigma^{4}}\left[\frac{y^{2}}{\sigma^{2}}-1\right] \exp \left(-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right) ∂y2∂2G(x,y)=2πσ41[σ2y2−1]exp(−2σ2x2+y2)

可得：
∇ 2 G ( x , y ) = − 1 π σ 4 [ 1 − x 2 + y 2 2 σ 2 ] exp ⁡ ( − x 2 + y 2 2 σ 2 ) \nabla^{2} G(x, y)=-\frac{1}{\pi \sigma^{4}}\left[1-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right] \exp \left(-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right) ∇2G(x,y)=−πσ41[1−2σ2x2+y2]exp(−2σ2x2+y2)

称为高斯-拉普拉斯算子（Laplacian of a Gaussian）简称LoG算子，也称为Marr边缘检测算子。

应用LoG算子时，高斯函数中标准差参数 σ σ σ的选择很关键，对图像边缘检测效果有很大的影响，对于不同图像应选择不同参数。

1.2.6 Canny [5]

Canny边缘检测的基本思想就是首先对图像选择一定的Gauss滤波器进行平滑滤波，然后采用非极值抑制技术进行处理得到最后的边缘图像。其步骤如下：

（1）用高斯滤波器平滑图像，减少噪声。
S ( x , y ) = g ( x , y ) ∗ G ( x , y ) S(x, y)=g(x, y)^{*} G(x, y) S(x,y)=g(x,y)∗G(x,y)

其中 g ( x , y ) g(x,y) g(x,y)为图像数据， G ( x , y ) G(x,y) G(x,y)为高斯函数。

（2）用一阶偏导的有限差分来计算梯度的幅值和方向。使用一阶差分卷积模板对图像进行卷积：
H 1 = ∣ − 1 − 1 1 1 ∣ , H 2 = ∣ 1 − 1 1 − 1 ∣ H_{1}=\left|\begin{array}{cc}-1 & -1 \\1 & 1\end{array}\right|, \quad H_{2}=\left|\begin{array}{cc}1 & -1 \\1 & -1\end{array}\right| H1=∣∣∣∣−11−11∣∣∣∣,H2=∣∣∣∣11−1−1∣∣∣∣
φ 1 ( x , y ) = S ( x , y ) ∗ H 1 ( x , y ) , φ 2 ( x , y ) = S ( x , y ) ∗ H 2 ( x , y ) \varphi_{1}(x, y)=S(x, y)^{*} H_{1}(x, y), \quad \varphi_{2}(x, y)=S(x, y)^{*} H_{2}(x, y) φ1(x,y)=S(x,y)∗H1(x,y),φ2(x,y)=S(x,y)∗H2(x,y)

则可得到每个点的局部梯度：
φ ( x , y ) = φ 1 2 ( x , y ) + φ 2 2 ( x , y ) \varphi(x, y)=\sqrt{\varphi_{1}^{2}(x, y)+\varphi_{2}^{2}(x, y)} φ(x,y)=φ12(x,y)+φ22(x,y)

和梯度方向
θ = tan ⁡ − 1 φ 2 ( x , y ) φ 1 ( x , y ) \theta=\tan ^{-1} \frac{\varphi_{2}(x, y)}{\varphi_{1}(x, y)} θ=tan−1φ1(x,y)φ2(x,y)
梯度方向与边缘方向关系：

（3）对梯度幅值进行非极大值抑制得到细化的边缘。

步骤（2）中确定的边缘点产生梯度边缘方向中的脊线。算法沿这些脊线（梯度方向）的顶部进行追踪，并将实际上不在脊线顶部的像素设置为零，从而在输出中给出一条细线，这个过程称为非极大值抑制。
如图所示，4个等角度的扇区的标号为0到3，对应3×3邻域的4种可能组合(比如：1-5,2-6,3-7,4-8)。

计算每个像素点的梯度方向 θ θ θ，根据 θ θ θ所在的扇区区间，让该像素与对应扇区区间的两个像素相比。如果该点的梯度值不比沿梯度线的两个相邻像素梯度值大，则令该像素为0。

（4）用双阈值算法检测和连接边缘。
使用两个阈值 T 1 T_1 T1和 T 2 T_2 T2 （ T 1 < T 2 T_1<T_2 T1<T2 ），其值大于 T 2 T_2 T2的脊线像素称为“强”边缘像素( B W 1 BW1 BW1)，值在 T 1 T_1 T1和 T 2 T_2 T2之间的脊线像素称为“弱”边缘像素( B W 2 BW2 BW2)。

由于 B W 1 BW1 BW1使用高阈值得到，因而含有很少的假边缘，但有间断（不闭合）。双阈值法要在 B W 1 BW1 BW1中把边缘连接成轮廓，当到达轮廓的端点时，该算法就在 B W 2 BW2 BW2的8邻点中寻找可以连接到轮廓上的边缘，试法将 B W 1 BW1 BW1的间断连接起来为止。
也就是说， T 2 T_2 T2用来找到每条线段， T 1 T_1 T1用来连接间断。

2 上代码！

2.1 边界处理函数

function imgOut = padding(imgIn,fsize,type)
% 函数功能：边界处理函数
% 关注微信公众号【二进制人工智能】，回复byjc获取demo全部代码
% 输入：
%      imgIn:输入图片
%      fsize：模板大小
%      type:padding类型
% 输出：
%      imgOut：处理后的图像
for i=1:(fsize-1)/2  % 模板为3，填充一层 模板为5填充两层，以此类推[m,n]=size(imgIn);imgTem=zeros(m+2,n+2);switch(lower(type))case 'zeros'        imgTem(2:m+1,2:n+1)=imgIn;case 'borden_value'  imgTem(2:m+1,2:n+1)=imgIn;% 边沿imgTem(1,2:n+1)=imgIn(1,:); imgTem(m+2,2:n+1)=imgIn(m,:); imgTem(2:m+1,1)=imgIn(:,1); imgTem(2:m+1,n+2)=imgIn(:,n);% 边角imgTem(1,1)=imgIn(1,1); imgTem(1,n+2)=imgIn(1,n);imgTem(m+2,1)=imgIn(m,1);imgTem(m+2,n+2)=imgIn(m,n);otherwisedisp('please input ''zeros''/''borden_value''') break;endimgIn=imgTem;
end
imgOut=imgIn;
end

2.2 Canny算法函数

function imgOut = cannyMyself(imgIn,thresholdValue,sigma)
% 函数功能：canny算法
% 关注微信公众号【二进制人工智能】，回复【byjc】获取demo全部代码
% 输入：
%      imgIn: 输入图片
%      thresholdValue：阈值
%      sigma: 高斯滤波标准差
% 输出：
%      imgOut：处理后的二值图像
imgIn = double(imgIn);
[m,n] = size(imgIn);
% 高斯滤波直接调用MATLAB中已有的函数
gaussFilter = fspecial('gaussian', [5,5], sigma);
imgIn = imfilter(imgIn,gaussFilter,'replicate');
wX=[-1,-1;1,1]/4;
wY=[1,-1;1,-1]/4;
% 梯度方向确定
theta = zeros(m, n);
sector = zeros(m, n);
imgOut=zeros(size(imgIn));
%标记：
% 2 1 0
% 3 X 3
% 0 1 2
for r = 1:(m-1)for c = 1:(n-1)childImg = imgIn(r:r+1,c:c+1);gx = wX .* childImg;gx = sum(gx(:));gy = wY .* childImg;gy = sum(gy(:));imgOut(r,c) = (gx^2+gy^2)^0.5 ;theta(r,c) = atand(gy/gx)  ;tem = theta(r,c);if (tem>22.5)&&(tem<67.5)sector(r,c) =  0;    elseif (tem<22.5)&&(tem>-22.5)sector(r,c) =  3;    elseif (tem>-67.5)&&(tem<-22.5)sector(r,c) = 2;    elsesector(r,c) = 1;    endend
end
%非极大值抑制
imgOut = [zeros(m,1),imgOut,zeros(m,1)];
imgOut = [zeros(1,n+2);imgOut;zeros(1,n+2)];
for r = 2:m+1for c = 2:n+1switch sector(r-1,c-1)case 0if((imgOut(r,c)<imgOut(r+1,c-1))| (imgOut(r,c)<imgOut(r-1,c+1)))bw1(r-1,c-1) = 0;elsebw1(r-1,c-1) = imgOut(r,c);endcase 1if((imgOut(r,c)<imgOut(r,c-1))| (imgOut(r,c)<imgOut(r,c+1)))bw1(r-1,c-1) = 0;elsebw1(r-1,c-1) = imgOut(r,c);endcase 2if((imgOut(r,c)<imgOut(r-1,c-1))| (imgOut(r,c)<imgOut(r+1,c+1)))bw1(r-1,c-1) = 0;elsebw1(r-1,c-1) = imgOut(r,c);endcase 3if((imgOut(r,c)<imgOut(r+1,c))| (imgOut(r,c)<imgOut(r-1,c)))bw1(r-1,c-1) = 0;elsebw1(r-1,c-1) = imgOut(r,c);endendend
end
for r = 1:mfor c = 1:n% 小阈值确定的边缘if (bw1(r,c)>thresholdValue(1))bw2(r,c) = 1;elsebw2(r,c) = 0;end% 大阈值确定的边缘if(bw1(r,c)>thresholdValue(2))bw3(r,c)=1;elsebw3(r,c)=0;endend
end
% bw2补充bw3
imgOut = bw3;
for r = 2:m-1for c = 2: n-1if (bw2(r,c)==1 & (imgOut(r,c)==1|imgOut(r+1,c)==1 ...|imgOut(r+1,c-1)==1|imgOut(r+1,c+1)==1|imgOut(r,c-1)==1 ...|imgOut(r,c+1)==1|imgOut(r-1,c-1)==1|imgOut(r-1,c)==1 ...|imgOut(r-1,c+1)==1))  % b2为边缘点且其所在邻域中b3也有边缘点imgOut(r,c) = 1;elseimgOut(r,c) = 0;endend
end

2.3 边缘检测函数

function imgOut=edgeMyself(imgIn,type,thresholdValue,sigma)
% 函数功能：边缘检测
% 关注微信公众号【二进制人工智能】，回复【byjc】获取demo全部代码
% 输入：
%      imgIn: 输入图片
%      type：边缘检测类型
%      thresholdValue：阈值
%      sigma: 高斯滤波标准差
% 输出：
%      imgOut：处理后的二值图像
if ndims(imgIn)==3imgIn = rgb2gray(imgIn);
end
imgIn=double(imgIn);
imgIn=imgIn/255;        % 像素值归一化到[0,1]
[m,n] = size(imgIn);
imgOut = imgIn;
if nargin<4sigma=1;
end
switch(lower(type))case 'roberts'wSize=2;wX = [-1,0;0,1]/2;wY = [0,-1;1,0]/2;case 'prewitt'wSize=3;wX=[1,1,1;0,0,0;-1,-1,-1]/6;wY=[1,0,-1;1,0,-1;1,0,-1]/6;imgIn=padding(imgIn,3,'borden_value'); case 'sobel'wSize=3;wX=[1,2,1;0,0,0;-1,-2,-1]/8;wY=[1,0,-1;2,0,-2;1,0,-1]/8;imgIn=padding(imgIn,3,'borden_value'); case 'laplacian'wSize=3;w=[-1,-1,-1;-1,8,-1;-1,-1,-1];imgIn=padding(imgIn,3,'borden_value'); case 'log'wSize=3;w=[-1,-1,-1;-1,8,-1;-1,-1,-1];% 高斯滤波直接调用MATLAB中已有的函数，具体实现可参考我的另一篇文章：%https://mp.weixin.qq.com/s?__biz=Mzg2MzM4Njk4MA==&mid=2247486698&idx=1&sn=9ac0d5c2e34b604bf04bbfa5b718b0f2&chksm=ce7820e1f90fa9f7830f7c6c222e9e032d63b86725d368415eea0f81f7322bff432f0246fd8b&token=2098489926&lang=zh_CN#rdwGauss = fspecial('gaussian',sigma);imgIn = imfilter(imgIn,wGauss,'replicate');imgIn=padding(imgIn,3,'borden_value'); case 'canny'imgOut=cannyMyself(imgIn,thresholdValue,sigma);returnotherwisedisp('please input roberts|prewitt|sobel|laplacian|log|canny')return
end
for r=1:mfor c=1:n% 模板大小为2时原图像无需padding处理，导致循环中需要特殊对待if wSize==2if c==n||r==mcontinueendchildImg = imgIn(r:r+1,c:c+1); elser_=r+1;c_=c+1;childImg = imgIn(r_-1:r_+1,c_-1:c_+1);endswitch(lower(type))case 'laplacian'   gradient=w.*childImg;gradient = abs(sum(gradient(:)));case 'log'gradient=w.*childImg;gradient = abs(sum(gradient(:)));otherwise             gradientX = wX .* childImg;gradientX = abs(sum(gradientX(:)));gradientY = wY .* childImg;gradientY = abs(sum(gradientY(:)));% 水平和垂直方向卷积结果之和作为梯度值gradient=gradientX+gradientY;endif(gradient > thresholdValue)imgOut(r,c) = 1;elseimgOut(r,c) = 0;end end
end
end

3 demo结果

原图：

边缘检测结果：

参考：

[1]https://blog.csdn.net/zaishuiyifangxym/article/details/89840396

[2]https://www.itdaan.com/blog/2016/06/29/c829f1a57ca24ee3a723bf236924aa85.html

[3]https://blog.csdn.net/aidem_brown/article/details/81542996

[4]https://blog.csdn.net/weixin_44378835/article/details/109406393

[5] https://blog.csdn.net/weixin_44378835/article/details/110572455

[6] https://blog.csdn.net/Prince_IT/article/details/90205259

[7] https://blog.csdn.net/humanking7/article/details/46606791