Canny算子作边缘检测

Canny边缘检测算子是John F. Canny于 1986 年开发出来的一个多级边缘检测算法。更为重要的是 Canny 创立了边缘检测计算理论(Computational theory of edge detection)解释这项技术如何工作。

文章结尾附源代码，可直接使用，有注释，便于理解。

算法发展

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

好的检测- 算法能够尽可能多地标识出图像中的实际边缘。

好的定位- 标识出的边缘要尽可能与实际图像中的实际边缘尽可能接近。

最小响应- 图像中的边缘只能标识一次，并且可能存在的图像噪声不应标识为边缘。

为了满足这些要求 Canny 使用了变分法，这是一种寻找满足特定功能的函数的方法。最优检测使用四个指数函数项的和表示，但是它非常近似于高斯函数的一阶导数。

算法步骤

去噪声

任何边缘检测算法都不可能在未经处理的原始数据上很好地处理，所以第一步是对原始数据与高斯平滑模板作卷积，得到的图像与原始图像相比有些轻微的模糊(blurred)。这样，单独的一个像素噪声在经过高斯平滑的图像上变得几乎没有影响。

寻找图像中的亮度梯度

图像中的边缘可能会指向不同的方向，所以 Canny 算法使用 4 个 mask 检测水平、垂直以及对角线方向的边缘。原始图像与每个 mask 所作的卷积都存储起来。对于每个点我们都标识在这个点上的最大值以及生成的边缘的方向。这样我们就从原始图像生成了图像中每个点亮度梯度图以及亮度梯度的方向。

在图像中跟踪边缘

较高的亮度梯度比较有可能是边缘，但是没有一个确切的值来限定多大的亮度梯度是边缘多大，所以 Canny 使用了滞后阈值。

滞后阈值需要两个阈值--高阈值与低阈值。假设图像中的重要边缘都是连续的曲线，这样我们就可以跟踪给定曲线中模糊的部分，并且避免将没有组成曲线的噪声像素当成边缘。所以我们从一个较大的阈值开始，这将标识出我们比较确信的真实边缘，使用前面导出的方向信息，我们从这些真正的边缘开始在图像中跟踪整个的边缘。在跟踪的时候，我们使用一个较小的阈值，这样就可以跟踪曲线的模糊部分直到我们回到起点。

一旦这个过程完成，我们就得到了一个二值图像，每点表示是否是一个边缘点。

一个获得亚像素精度边缘的改进实现是在梯度方向检测二阶方向导数的过零点。

它在梯度方向的三阶方向导数满足符号条件

...表示用高斯核平滑原始图像得到的尺度空间表示 L 计算得到的偏导数。用这种方法得到的边缘片断是连续曲线，这样就不需要另外的边缘跟踪改进。滞后阈值也可以用于亚像素边缘检测。

算法参数

Canny 算法包含许多可以调整的参数，它们将影响到算法的计算的时间与实效。

高斯滤波器的大小:第一步所用的平滑滤波器将会直接影响 Canny 算法的结果。较小的滤波器产生的模糊效果也较少，这样就可以检测较小、变化明显的细线。较大的滤波器产生的模糊效果也较多，将较大的一块图像区域涂成一个特定点的颜色值。这样带来的结果就是对于检测较大、平滑的边缘更加有用，例如彩虹的边缘。

阈值:使用两个阈值比使用一个阈值更加灵活，但是它还是有阈值存在的共性问题。设置的阈值过高，可能会漏掉重要信息;阈值过低，将会把枝节信息看得很重要。很难给出一个适用于所有图像的通用阈值。目前还没有一个经过验证的实现方法。

算法结论

Canny 算法适用于不同的场合。它的参数允许根据不同实现的特定要求进行调整以识别不同的边缘特性。对于PC上的实时图像处理来说可能慢得无法使用，尤其是在使用大的高斯滤波器的情况下。但是，我们讨论计算能力的时候，也要考虑到随着处理器速度不断提升，有望在未来几年使得这不再成为一个问题。

#图像数据可视化模块
import matplotlib.pyplot as plt
#图像色彩映射模块
import matplotlib.cm as cm
#算子与图像矩阵处理模块
import numpy as np
#数学计算和数学常量处理模块
import math

#载入原图
img = plt.imread(r"C:\Users\86159\PycharmProjects\pythonProject2\test.jpg.jpg")
#设置高斯滤波器标准差，缺省值为1
sigma1 = sigma2 =1
sum = 0

#初始化5*5的高斯算子矩阵
gaussian = np.zeros([5,5])
for i in range(5):
    for j in range(5):
        #生存二维高斯分布矩阵
        gaussian[i,j]=math.exp(-1/2*(np.square(i-3/np.square(sigma1)+np.square(j-3)/np.square(sigma2))))/(2*math.pi*sigma1*sigma1)
        sum = sum + gaussian[i,j]
gaussian = gaussian/sum

def rgb2gray(rgb):
    #图像灰度化处理
    return np.dot(rgb[...,:3],[0.299,0.587,0.114])

#高斯滤波
gray = rgb2gray(img)

W , H =gray.shape
new_gray = np.zeros([W-5,H-5])

for i in range(W-5):
    for j in range(H-5):
        #与高斯矩阵卷积实现滤波
        new_gray[i,j]=np.sum(gray[i:i+5,j:j+5]*gaussian)

#通过求梯度幅值使图像增强
W1,H1=new_gray.shape
dx = np.zeros([W1-1,H1-1])
dy = np.zeros([W1-1,H1-1])
d  =np.zeros([W1-1,H1-1])

for i in range(W1-1):
    for j in range(H1-1):
        dx[i,j]=new_gray[i,j+1]-new_gray[i,j]
        dy[i,j]=new_gray[i+1,j]-new_gray[i,j]
        #图像梯度幅值作为图像增强值
        d[i,j]=np.sqrt(np.square(dx[i,j])+np.square(dy[i,j]))

#非极大值抑制NMS
W2,H2 = d.shape
NMS = np.copy(d)
NMS[0,:]=NMS[W2-1,:]=NMS[:,0]=NMS[:,H2-1]=0

for i in range(1,W2-1):
    for j in range(1,H2-1):
        if d[i,j] == 0:
            NMS[i,j] = 0
        else:
            gradX = dx[i,j]
            gradY = dy[i,j]
            gradtemp = d[i,j]

            #Y方向梯度值较大：
            if np.abs(gradY) > np.abs(gradX):
                weight = np.abs(gradX)/np.abs(gradY)
                grad2 = d[i-1,j]
                grad4 = d[i+1,j]

                #方向梯度值符号相同
                if gradX * gradY >0:
                    grad1 = d[i-1,j-1]
                    grad3 = d[i+1,j+1]
                #方向梯度值符号相反
                else:
                    grad1 = d[i-1,j+1]
                    grad3 = d[i+1,j-1]
            #X方向梯度值较大
            else:
                weight = np.abs(gradY)/np.abs(gradX)
                grad2=d[i,j-1]
                grad4=d[i,j+1]

                if gradX * gradY >0:
                    grad1 = d[i+1,j-1]
                    grad3 = d[i-1,j+1]
                else:
                    grad1 = d[i-1,j-1]
                    grad3 = d[i+1,j+1]

                gradtemp1 = weight*grad1+(1-weight)*grad2
                gradtemp2 = weight*grad3+(1-weight)*grad4
                if gradtemp >= gradtemp1 and gradtemp>=gradtemp2:
                    NMS[i,j]=gradtemp
                else:
                    NMS[i,j]=0

#双阈值算法检测，连接边缘
W3,H3 = NMS.shape
DT = np.zeros([W3,H3])

#低阈值TL，高阈值TH
TL = 0.2*np.max(NMS)
TH = 0.3*np.max(NMS)

for i in range(1,W3-1):
    for j in range(1,H3-1):
        if (NMS[i,j]<TL):
            DT[i,j]=0
        elif (NMS[i,j]>TH):
            DT[i,j]=1
        elif (NMS[i-1,j-1:j+1]<TH).any() or (NMS[i+1,j-1:j+1]).any() or (NMS[i,[j-1,j+1]]<TH).any():
            DT[i,j]=1

plt.subplot(2,2,1)
#从原图像转化的灰度图像
plt.imshow(new_gray,cmap=cm.gray)
plt.axis('off')
plt.subplot(2,2,2)
#高斯滤波后的图像
plt.imshow(d,cmap=cm.gray)
plt.axis('off')
plt.subplot(2,2,3)
#非极大值抑制图像
plt.imshow(NMS,cmap=cm.gray)
plt.axis('off')
plt.subplot(2,2,4)
#双阈值检测边缘图像
plt.imshow(DT,cmap=cm.gray)
plt.axis('off')
plt.show()