CV06_Canny边缘检测算法和python实现

1.1简介

Canny边缘检测算法是计算机视觉和图像处理领域中一种广泛应用的边缘检测技术，由约翰·F·坎尼（John F. Canny）于1986年提出。它是基于多级处理的边缘检测方法，旨在实现以下三个优化目标：

好的检测：尽可能多地检测出真正的边缘，同时尽量减少假阳性（误报）。
好的定位：检测到的边缘应该尽可能接近真实边缘的实际位置。
最小响应：对于单个边缘，只希望有一个响应，避免重复检测同一边缘。

Canny算法的步骤可以概括为以下几个阶段：

1. 高斯滤波（降噪）

算法的第一步是对原始图像进行高斯滤波，以去除图像中的噪声。高斯滤波是一种平滑处理，通过与高斯核进行卷积来实现，可以有效减少图像中的随机波动，而不模糊边缘太多。

2. 计算梯度幅度和方向

接着，算法会计算图像中每个像素点的梯度幅度和方向。梯度代表了像素值变化的强度，是边缘检测的关键。这通常通过应用Sobel算子或类似的微分算子在水平和垂直方向上进行，然后组合这些结果来得到梯度的大小和方向。梯度方向用于后续的非极大值抑制步骤。

3. 非极大值抑制（Edge thinning）

在这一步，算法沿着每个像素的梯度方向检查周围的像素，仅保留局部梯度最大的像素，即真正的边缘候选点。这样做是为了细化边缘，确保每个边缘只有单一的一行像素宽度，同时进一步减少噪声。

4. 双阈值检测和边缘连接

Canny算法采用了滞后阈值策略来区分强边缘和弱边缘。首先设定两个阈值（高阈值和低阈值），梯度幅度大于高阈值的像素被标记为强边缘，小于低阈值的像素被视为背景。介于两者之间的像素，如果与已标记的强边缘像素相连，则被视为弱边缘并最终加入到边缘集合中。这样可以有效连接断裂的边缘，同时排除孤立的噪声点。

5. 边缘跟踪与抑制

最后，通过边缘跟踪和抑制进一步优化边缘，确保边缘连续且光滑，同时减少不必要的边缘响应。

Canny算法之所以被广泛认为是最优的边缘检测算法之一，是因为它在准确性、效率和鲁棒性之间取得了很好的平衡。尽管如此，它也并非没有缺点，例如计算成本相对较高，且对阈值的选择敏感。但通过适当的参数调整和优化，Canny算法在许多应用场景中仍能提供高质量的边缘检测效果。

1.2 差分和梯度的关系

离散值的梯度计算通常涉及到差分，特别是在数字图像处理、数值分析以及某些类型的机器学习算法中。梯度本质上是一个向量，表示函数在某一点上的方向导数，即函数增长最快的方向。在连续空间中，梯度是偏导数的集合。但在离散空间（如像素网格上的图像）中，我们不能直接计算导数，而是通过有限差分近似来估计梯度。

一维情况

对于一维离散函数 𝑓(𝑥)，在点 𝑥𝑖处的梯度（即导数的近似）可以通过向前差分或向后差分来估算：

这里，ℎ 是两点之间的间隔，通常代表采样间隔。

二维情况

在二维图像处理中，一个像素点 (𝑥𝑖,𝑦𝑗) 处的梯度是一个二维向量 (∂𝑓∂𝑥,∂𝑓∂𝑦)，分别表示函数在 𝑥 和 𝑦 方向上变化的速率。同样使用差分近似：

最终，梯度的模（即梯度的大小）和方向常被用于边缘检测、图像增强、特征提取等应用中。在某些情况下，为了得到连续且更平滑的梯度估计，可以使用更复杂的方法，如高斯差分或其他平滑先验的梯度计算方法。

为什么离散函数的梯度用差分来表示

离散函数的梯度计算之所以使用差分来表示，是因为在离散空间中（比如数字图像中的像素点），我们无法直接应用连续函数的微分概念。在连续空间中，导数定义为函数值的变化率，即无限小的区间内的变化量比上这个区间的长度趋向于零的极限。然而，在离散空间里，函数值只在离散的点上有定义，没有“无限小”的概念，因此无法直接求极限。

差分则是离散空间中导数概念的自然近似。它通过计算相邻两点间函数值之差来量化函数值的变化率。这与导数的思想相似，但适用于离散数据点。例如，在一维离散函数中，两点之间的差分可以看作是在这两个点之间函数斜率的近似。在二维图像处理中，通过计算像素间灰度值的差异，我们可以近似得到图像中亮度变化的速率和方向，从而识别边缘、纹理等特征。

因此，差分不仅提供了计算离散梯度的有效方法，而且在数值分析和图像处理等领域的实际应用中非常实用，能够以一种直观且计算上可行的方式处理和分析离散数据的局部变化特性。

1.3 三种邻域像素距离定义

Canny算法中使用欧氏距离和城市距离都可以，欧氏距离含根号计算比较慢但是比较精准，城市距离都是整数计算比较快。

1.4 Sobel(索贝尔)算子

索贝尔算子（Sobel Operator）是一种广泛应用于图像处理领域的边缘检测技术，特别适用于需要确定图像中物体边界的位置和方向的情况。它是基于离散微分的一种方法，通过对图像进行局部微分操作，来估计图像中每个像素点的梯度强度及其方向，进而发现边缘。

原理与特点

离散差分算子：索贝尔算子是一种一阶离散差分算子，能够提供图像亮度函数梯度的近似值。它通过计算像素点在水平和垂直方向上的梯度分量来工作，这两个分量分别是图像在x轴和y轴方向上的变化率。
权重分配：与简单的梯度估计方法相比，索贝尔算子在计算像素差分时考虑了邻域像素的加权平均，即距离当前像素越近的邻像素被赋予的权重越大。这种权重分配机制使得索贝尔算子对噪声具有一定的抑制作用，因为它不仅仅依赖于直接相邻的像素差异，而是综合考虑了一个小邻域内的变化。
模板结构：索贝尔算子通常使用两个3x3的卷积核（掩模）来实现，一个用于计算水平梯度（Sobel x方向），另一个用于计算垂直梯度（Sobel y方向）。这些模板的系数经过精心设计，既能捕捉到方向上的变化，又能保持对噪声的一定鲁棒性。
边缘检测与方向：通过计算每个像素点的梯度幅度（通常是x和y方向梯度平方和的平方根），可以得到边缘强度。同时，梯度的方向可以指示边缘的方向。
性能与局限：索贝尔算子因其易于实现、计算效率高且能提供较好的边缘定位而被广泛采用。但它也存在一些局限性，比如检测到的边缘可能会比较粗，并且在某些情况下可能产生伪边缘，特别是在图像噪声较大或者边缘不清晰的情况下。

1.3 算法思路

降噪（滤波）

在对图像进行边缘检测之前我们需要对图像进行滤波，常用的一般是高斯滤波，这样能够过滤点一些噪声点，否则在后期边缘检测会出现一些斑点。

滤波也有缺点，它会使图像变得模糊，从而丢失掉一些细节信息，也会导致边缘线一定程度的变粗。

灰度化图像并计算梯度

首先，我们的RGB图需要转成灰度图才能进行边缘检测，灰度图每个点有一个灰度值，最亮是255，最暗是0。边缘和里外两侧灰度值变化比较陡峭，但边缘连线上的灰度值变化比较平滑，所以边缘都是一些突变点。

我们用亮度值来衡量边缘，把变化越剧烈的点变得越亮，变化越平缓的点越暗，最终就能得到边缘检测的效果。

对于一维的离散函数，我们用梯度差分来表示亮度值，即 Dx = Ax+1 - Ax。比如我们有一列离散值1，2，10，100，100，这些是灰度值。对它们求差分，得到1，8，90，0，这些是亮度值。显然亮度值为90的点是边缘。

对于二维的离散函数，就比如图像的像素值，我们可以用前向水平和垂直差分表示：

$\overrightarrow{D}=\begin{bmatrix}D^{\prime}x\\D^{\prime}y\end{bmatrix}=\begin{bmatrix}a_{x+1,y}-a_{x,y}\\a_{x,y+1}-a_{x,y}\end{bmatrix}$

D的模值就当作亮度值，写作 $\overrightarrow{D}=\sqrt{D_{x}^{2}+D_{y}^{2}}$ ，D的方向角θ，记作 $\theta=\arctan(\frac{D'y}{D'x})$

但是我们在Canny边缘检测时，梯度用的不是前向差分而是中心差分：

$\overrightarrow{D}_{xy}=\begin{bmatrix}D_{'x}\\D_{'y}\end{bmatrix}=\begin{bmatrix}\frac{a_{x+1,y}-a_{x-1,y}}{2}\\\\\frac{a_{x,y+1}-a_{x,y-1}}{2}\end{bmatrix}$

a1	a2	a3
a4	a5	a6
a7	a8	a9

但是我们可以发现，用这种方法算中心点a5的梯度，a1,a3,a7,a9是没有作用的，但实际上它们对边缘检测是有影响的。这就引入了Sobel算子。

对于a5使用索贝尔算子后，采用城市距离，的中心差分梯度为：

$\left.Da_{5}=\left[\begin{matrix}D_{5}^{\prime}x\\D_{5}^{\prime}y\end{matrix}\right.\right]=\left[\begin{matrix}\frac{a_{6}-a_{4}}{2\times1}+\frac{a_{9}-a_{1}}{2\times2}&-\frac{a_{7}-a_{3}}{2\times2}\\\frac{a_{8}-a_{2}}{2\times1}&+\frac{a_{9}-a_{1}}{2\times2}&+\frac{a_{7}-a_{3}}{2\times2}\end{matrix}\right]$ ，计算求模再向下取整。

公式如下：

这里第一行最后一个是减法的原因是根据sobel算子和XY轴正方向，a7到a3的方向在X轴的分量是负的。

得到θ角（范围-180°，180°）以后，我们规定X轴正方向为0°。假如θ角是45度，那么说明从a5到a9这个方向亮度值变化最剧烈。

计算图像中每一个像素点的亮度值（梯度），比较像素点与梯度方向相邻点的梯度值，最大的那个即为边缘。

非极大值抑制

如果只求梯度会导致一个问题，就是导致求出的图像边缘比较模糊而且线条比较暗且线条比较粗。

我们放大原图可以看到，边缘有几行像素只有最中间是最黑的，我们假设灰度值为0，旁边的像素灰度值可能是10，15，20，也就是说，边缘到非边缘其实是一个渐进的变化，不是突然从0到200的，这就导致了边缘的线条比较粗。

如何解决这个问题？就是用NMS。我们再看这个图，我们已经找到了最大梯度值的边缘，但是其他的值得到了保留，例如第三列的7，上面的5和下面的4还是保留着的，这时候我们就用NMS，将他们置为零，这样真正的边缘就出现了。

可以明显看到线条变细了。

双阈值检测

比如上面的图，如果某些像素点灰度值大于180,那么标记为强边缘，直接置为255，相对的，比如某些像素点的灰度值小于70，标记为背景，那么直接置为0.

孤立弱边缘抑制

什么是弱边缘，就是灰度值介于背景和强边缘之间的那些像素点的集合。弱边缘分为有用的和无用的，有用的大多与强边缘相近，无用的大多是孤立的。所以，我们要抑制这种孤立的弱边缘。

那么怎么才算孤立呢？假如一个弱边缘像素点，它周围一圈八个点不含强边缘的像素点，那么这个像素点就是孤立的，我们就置为0.

思考

在计算sobel算子下的离散函数梯度时，图像最边缘的像素点缺少相邻像素，该怎么办？

在使用Sobel算子计算图像边缘的离散函数梯度时，确实会遇到边界问题，因为图像的边缘像素缺少完整的邻域来应用3x3的卷积核。处理这种情况有几种常见方法：

边界扩展（Padding）：
- 镜像扩展（Mirror Padding）：缺失的像素值由其对称位置的像素值镜像复制而来。
- 常数扩展（Constant Padding）：缺失的像素值填充为一个常数值，如0或图像边缘的像素值。
- 重复最后一个像素值（Replicate Padding）：边界像素值直接复制作为扩展。
- 反射扩展（Reflect Padding）：类似于镜像，但更像是物理反射，比如最后一个像素值的“前面”像素就是它自己。
忽略边界像素：
- 直接不计算图像边缘像素的梯度，这样会导致最终处理后的图像比原图小2个像素宽和高。
循环移位（Cyclic Padding）：
- 将图像视为周期性的，即图像的右侧与左侧相连，上边与下边相连，从而形成一个循环的边界条件。
使用专门针对边缘的算法变体：
- 有些算法对边界进行了特殊处理，设计了特殊的边界滤波器或者调整了边界处的卷积操作，以减少边缘效应。

在实际应用中，选择哪种方法取决于具体需求和对边缘效应的容忍度。镜像扩展是常用且效果较好的方法之一，因为它能够保持图像边缘的连续性，减少因边界处理引入的人工痕迹。在一些图像处理库中，如OpenCV，提供了对边界处理的选择项，允许用户根据具体情况选择合适的边界处理方式。

Python实现

1，rgb转灰度图

2，滤波（要注意滤波强度，滤波会使图像变模糊，图像太模糊会失去细节信息。但是滤波可以去除噪声点，也可以用边缘强化滤波强化边缘）

3，求梯度（sobel算子）

4，非极大值抑制（使边缘更清晰）

5，双阈值(消去噪点，强化边缘)

6，抑制孤立弱边缘

from PIL import Image
from matplotlib import pyplot as plt
import math
from PIL import ImageFilter

im = Image.open("pic.jpg")
im_gray = im.convert('L')   #rgb转灰度图
# im_gray.show()  #显示灰度图

# im_gray = im_gray.filter(ImageFilter.SMOOTH)

#Image转为list类型
im_array = im_gray.load()
im_list = [[0 for i in range(im.size[0])] for j in range(im.size[1])]
for i in range(im.size[1]):
    for j in range(im.size[0]):
        im_list[i][j] = im_array[j, i]
plt.figure(1)
plt.imshow(im_list, cmap='gray')
# plt.show()


#求梯度 （梯度模值和方向）
#方向：  0：0度(x)    1: 45度   2： 90度（y） 3: 135度
grad = [[[0, 0] for i in range(im.size[0])] for j in range(im.size[1])]  #170 * 200 * 2 每个像素点的梯度大小和方向
for j in range(1, im.size[1] - 1):
    for i in range(1, im.size[0] - 1):
        #x方向梯度
        grad_x = im_list[j + 1][i + 1] + im_list[j - 1][i + 1] + 2 * im_list[j][i + 1] - \
                 im_list[j - 1][i - 1] - im_list[j + 1][i - 1] - 2 * im_list[j][i - 1]
        #y方向梯度
        grad_y = im_list[j + 1][i - 1] + im_list[j + 1][i + 1] + 2 * im_list[j + 1][i] - \
                 im_list[j - 1][i - 1] - im_list[j - 1][i + 1] - 2 * im_list[j - 1][i]
        grad_x = math.floor(grad_x / 4)
        grad_y = math.floor(grad_y / 4)
        #合梯度
        grad[j][i][0] = math.floor(math.sqrt(grad_x * grad_x + grad_y * grad_y))
        if(grad[j][i][0] > 255):
            grad[j][i][0] = 255
        # if(grad[j][i][0] < 50):
        #     grad[j][i][0] = 0
        #梯度方向
        if(grad_x == 0):
            grad[j][i][1] = 2  #y方向
        else:
            theta = math.atan2(grad_y, grad_x)
            if(math.fabs(theta) < math.pi / 8):
                grad[j][i][1] = 0 #x方向
            elif(theta > 0):
                if(math.fabs(theta) <math.pi * 3 / 8):
                    grad[j][i][1] = 1  #45度方向
                else:
                    grad[j][i][1] = 2#y方向
            else:
                if (math.fabs(theta) < math.pi * 3 / 8):
                    grad[j][i][1] = 3  # 135度方向
                else:
                    grad[j][i][1] = 2  # y方向

#显示梯度图
img = [[0 for i in range(im.size[0])] for j in range(im.size[1])]
for i in range(im.size[1]):
    for j in range(im.size[0]):
        img[i][j] = grad[i][j][0]
plt.figure(2)
plt.imshow(img, cmap='gray')
# plt.show()


#非极大值抑制,使边缘更清晰和细
img2 = img
for j in range(1, im.size[1] - 1):
    for i in range(1, im.size[0] - 1):
        dir = grad[j][i][1]
        grad_now = grad[j][i][0]
        if (dir == 0):  #梯度方向为x
            if(grad_now < grad[j][i + 1][0] or grad_now < grad[j][i - 1][0]):
                # grad[j][i][0] == 0
                img2[j][i] = 0
                print('0')
        elif(dir == 1):#45度方向
            if(grad_now < grad[j + 1][i + 1][0] or grad_now < grad[j - 1][i - 1][0]):
                # grad[j][i][0] == 0
                img2[j][i] = 0
                print('1')
        elif(dir == 2): #y方向
            if(grad_now < grad[j + 1][i][0] or grad_now < grad[j - 1][i][0]):
                # grad[j][i][0] == 0
                img2[j][i] = 0
                print('2')
        else:  #145度
            if(grad_now < grad[j + 1][i - 1][0] or grad_now < grad[j - 1][i + 1][0]):
                # grad[j][i][0] == 0
                img2[j][i] = 0
                print('3')

#显示非极大值抑制后的图像
plt.figure(3)
plt.imshow(img2, cmap='gray')
# plt.show()




#双阈值检测
low = 30
high = 60

for i in range(im.size[1]):
    for j in range(im.size[0]):
        if(img2[i][j] < low):
            img2[i][j] = 0
        elif(img2[i][j] > high):
            img2[i][j] = 255

#双阈值检测后的图像
plt.figure(4)
plt.imshow(img2, cmap='gray')
# plt.show()



#抑制孤立的弱边缘
listx = [-1, 0, 1]
listy = [-1, 0, 1]
img3 = img2
for j in range(1, im.size[1] - 1):
    for i in range(1, im.size[0] - 1):
        flag = 1
        for dx in range(len(listx)):
            for dy in range(len(listy)):
                j = j + listy[dy]
                i = i + listx[dx]
                if(img2[j][i] == 255):
                    flag = 0
                    break
            if(not flag):
                break
        if(flag):
            img3[j][i] = 0

#抑制孤立的弱边缘后的图像
# print('ssssssss')
plt.figure(5)
plt.imshow(img3, cmap='gray')
plt.show()

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

(1)最优检测：算法能够尽可能多地标识出图像中的实际边缘，漏检真实边缘的概率和误检非边缘的概率都尽可能小；

(2)最优定位准则：检测到的边缘点的位置距离实际边缘点的位置最近，或者是由于噪声影响引起检测出的边缘偏离物体的真实边缘的程度最小；

(3)检测点与边缘点一一对应：算子检测的边缘点与实际边缘点应该是一一对应。