Python代码实现Canny算法——图像边缘轮廓提取

一、Canny算法主要思路：

1.通过高斯核进行滤波，去除噪声点，使图像平滑

2.通过求取偏导，使用梯度强度和方向来提取边缘轮廓

3.采用非极大值抑制的方法，将提取到的粗犷的边缘轮廓细化

4.采用双阈值来进一步过滤噪声即过滤掉不想要的线条，同时保证轮廓的连续性

1.高斯核滤波

众所周知，正态分布又叫做高斯分布，高斯核滤波其实就是正态分布的一个应用。

高斯分布：(和标准的二维正态分布有点差别)
$$G_{\sigma }=1/(2\pi {\sigma }_1{\sigma }_2)\ast e^{-(x-\mu {)}^2/{\sigma }_1^2+(y-\mu {)}^2)/{\sigma }_2^2}$$
我们假设一个3*3的高斯核：（我用Excel画的）

设高斯核中小格子的坐标为(x,y)， 高斯核中间的小格子为（0,0），左上角的小格子为（-1，-1），以此类推。

将高斯核中每个小格子的坐标带入到高斯分布公式中，就可以得到每个小格子的值。

高斯核中的所有值得核应该为1，因此需要进行归一化处理，即求出每个小格子占所有小格子总和的权重。

高斯核的应用——卷积

图像是以像素值来表示的。

如图为一个36*36的灰度图像，每个小格子中都有一个像素值，将高斯核覆盖到图像上，如图红色和黑色格子。计算图像和高斯核对应格子的值，并求其总和，即为黑色格子进行高斯卷积之后所对应的像素值。理论上计算时，通常会对高斯核进行上下镜像，之后再进行计算，但实际应用中通常不考虑该镜像操作。

以此类推，可计算整个图像中所有格子进行卷积后的值。再将这些值输出，即可得到进行高斯卷积之后的图像。

高斯卷积会使得图像丢失最外围的像素。因此通常会在卷积之前，在图像的周围填充一些像素，来保证图像卷积后大小不变。

2.偏导和梯度

为什么要求偏导或者梯度？

图像的边缘是由于像素值的突变引起的，因此我们需要通过求偏导或者梯度来获取到突变值发生的位置，提取轮廓。

（公式不好打，上图吧)

求偏导或者梯度函数，其本质上还是和高斯卷积是类似的。

常用的几个求偏导和梯度的模板：

Prewitt模板算子：
$$M_x=\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]M_y=\left[\begin{array}{ccc}1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\end{array}\right]$$
Sobel模板算子：
$$M_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]M_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]$$

3.非极大值抑制法

由于图像的像素是逐渐进行过渡的，因此我们通过上述过程提取到的边缘会很粗，我们想要得到的图像的边缘应该只有一个像素的粗细。因此我们要提取突变最大的地方作为图像的边缘。

非极大值抑制法，即沿梯度方向，选取梯度强度最大的一个像素块进行保留，将其余像素块舍弃的方法。

比较方法通常为：沿梯度方向，取一像素点，和该像素点前后各一像素的像素点，进行比较，保留最大值。

如图，假设该3*3的核的每一个格子均在所提出的轮廓线上。

假设梯度方向是沿x轴方向的，我们需要比较点（-1，0）和点（0，0）两处的梯度强度，若点（0，0）的梯度强度大于（-1，0），则比较（0,0）和（1,0）的梯度强度，若（0,0）小于（1,0），则继续比较（1,0）和（2,0），若（0,0）大于（-1,0）且大于（1,0）则（0,0）最大，可保留，其余的均丢弃，即像素值设为0.

然而实际应用中，梯度方向和边垂直，并不都是水平或者垂直的，因此我们可以选取其所靠近的值，采用不同的权重来计算该点的梯度强度。

如下图。

4.双阈值滤波

双阈值滤波是基于如下的假设进行的：

假设，所有噪声所产生的线条与图像的轮廓是不发生连接的。

这样，我们就可以设置两个阈值，一个高阈值，一个低阈值。低阈值用来过滤噪声和突变不明显的线条，高阈值用来保留轮廓线。在低阈值和高阈值之间，且和轮廓线相邻的线条可保留，否则舍弃。

二、代码实现

（代码仅供学习参考使用，尚有不足和待优化之处）

import math
import numpy as np
import matplotlib.pyplot as plt


# 生成高斯核
def gaussian_create():
    sigma1 = sigma2 = 1
    gaussian_sum = 0
    g = np.zeros([3, 3])
    for i in range(3):
        for j in range(3):
            g[i, j] = math.exp(-1 / 2 * (np.square(i - 1) / np.square(sigma1)
                                         + (np.square(j - 1) / np.square(sigma2)))) / (
                              2 * math.pi * sigma1 * sigma2)
            gaussian_sum = gaussian_sum + g[i, j]
    g = g / gaussian_sum  # 归一化
    return g


# 产生灰度图
def gray_fuc(rgb):
    return np.dot(rgb[..., :3], [0.299, 0.587, 0.114])


# 高斯卷积
def gaussian_blur(gray_img, g):
	'''
	gray_img:灰度图
	g:高斯核
	'''
    gray_img = np.pad(gray_img, ((1, 1), (1, 1)), constant_values=0)  # 填充
    h, w = gray_img.shape
    new_gray_img = np.zeros([h - 2, w - 2])
    for i in range(h - 2):
        for j in range(w - 2):
            new_gray_img[i, j] = np.sum(gray_img[i:i + 3, j:j + 3] * g)
    return new_gray_img


# 求高斯偏导
def partial_derivative(new_gray_img):
	'''
	new_gray_img:高斯卷积后的灰度图
	'''
    new_gray_img = np.pad(new_gray_img, ((0, 1), (0, 1)), constant_values=0)  # 填充
    h, w = new_gray_img.shape
    dx_gray = np.zeros([h - 1, w - 1])  # 用来存储x方向偏导
    dy_gray = np.zeros([h - 1, w - 1])  # 用来存储y方向偏导
    df_gray = np.zeros([h - 1, w - 1])  # 用来存储梯度强度
    for i in range(h - 1):
        for j in range(w - 1):
            dx_gray[i, j] = new_gray_img[i, j + 1] - new_gray_img[i, j]
            dy_gray[i, j] = new_gray_img[i + 1, j] - new_gray_img[i, j]
            df_gray[i, j] = np.sqrt(np.square(dx_gray[i, j]) + np.square(dy_gray[i, j]))
    return dx_gray, dy_gray, df_gray


# 非极大值抑制
def non_maximum_suppression(dx_gray, dy_gray, df_gray):
	'''
	dx_gray:x方向梯度矩阵
	dy_gray:y方向梯度矩阵
	df_gray:梯度强度矩阵
	'''
    df_gray = np.pad(df_gray, ((1, 1), (1, 1)), constant_values=0)  # 填充
    h, w = df_gray.shape
    for i in range(1, h - 1):
        for j in range(1, w - 1):
            if df_gray[i, j] != 0:
                gx = math.fabs(dx_gray[i - 1, j - 1])
                gy = math.fabs(dy_gray[i - 1, j - 1])
                if gx > gy:
                    weight = gy / gx
                    grad1 = df_gray[i + 1, j]
                    grad2 = df_gray[i - 1, j]
                    if gx * gy > 0:
                        grad3 = df_gray[i + 1, j + 1]
                        grad4 = df_gray[i - 1, j - 1]
                    else:
                        grad3 = df_gray[i + 1, j - 1]
                        grad4 = df_gray[i - 1, j + 1]
                else:
                    weight = gx / gy
                    grad1 = df_gray[i, j + 1]
                    grad2 = df_gray[i, j - 1]
                    if gx * gy > 0:
                        grad3 = df_gray[i + 1, j + 1]
                        grad4 = df_gray[i - 1, j - 1]
                    else:
                        grad3 = df_gray[i + 1, j - 1]
                        grad4 = df_gray[i - 1, j + 1]
                t1 = weight * grad1 + (1 - weight) * grad3
                t2 = weight * grad2 + (1 - weight) * grad4
                if df_gray[i, j] > t1 and df_gray[i, j] > t2:
                    df_gray[i, j] = df_gray[i, j]
                else:
                    df_gray[i, j] = 0
    return df_gray


# 双阈值过滤
def double_threshold(df_gray, low, high):
	'''
	df_gray:梯度强度矩阵
	low:低阈值
	high:高阈值
	'''
    h, w = df_gray.shape
    for i in range(1, h - 1):
        for j in range(1, w - 1):
            if df_gray[i, j] < low:
                df_gray[i, j] = 0
            elif df_gray[i, j] > high:
                df_gray[i, j] = 1
            elif (df_gray[i, j - 1] > high) or (df_gray[i - 1, j - 1] > high) or (
                    df_gray[i + 1, j - 1] > high) or (df_gray[i - 1, j] > high) or (df_gray[i + 1, j] > high) or (
                    df_gray[i - 1, j + 1] > high) or (df_gray[i, j + 1] > high) or (df_gray[i + 1, j + 1] > high):
                df_gray[i, j] = 1
            else:
                df_gray[i, j] = 0
    return df_gray


if __name__ == '__main__':
    # 读取图像
    img = plt.imread('C:/Users/Administrator/Desktop/zhongqiu.jpg')
    # 生成高斯核
    gaussian = gaussian_create()
    # 生成灰度图
    gray = gray_fuc(img)
    # 高斯卷积
    new_gray = gaussian_blur(gray, gaussian)
    # 求偏导
    d = partial_derivative(new_gray)
    dx = d[0]
    dy = d[1]
    df = d[2]
    # 非极大值抑制
    new_df = non_maximum_suppression(dx, dy, df)
    # 双阈值过滤,并将图像转换成转化二值图
    low_threshold = 0.15 * np.max(new_df)
    high_threshold = 0.2 * np.max(new_df)
    result = double_threshold(new_df, low_threshold, high_threshold)
    # 输出图像
    plt.imshow(result, cmap="gray")
    plt.axis("off")
    plt.show()

运行结果

（不晓得为啥上次上传的图片显示违规。。。换了张图）上图为原图，下图为提取到的轮廓。运行时间大约44s

更换图片，则要根据需要调整双阈值滤波的高低阈值大小，来获得满意的输出。

代码改进

上述代码中采用的双阈值过滤函数，是判断弱边周围是否存在强边，来进而确定该弱边是否是我们所需要的边，是否进行滤除。该方法的缺陷：

强边为H，其上有一点a与弱边相连。该弱边为L，其上有一点b和一点c，b点和a点相连。
如果，在判断弱边是否滤除的时候，先判断b点，后判断c点，得知b点和a点相连，b点设为1值保留，c点和b点相连，因此c点也设为1值保留
如果，先判断c点，后判断b点，那么得出的结论是：c点周围没有1值，丢弃。这将导致我们丢失掉我们想要的点

因此改变双阈值过滤函数的算法思想，判断强边周围是否存在弱边，即通过强边延伸弱边。这样可以使我们提取到的边更加的完整。

# 双阈值过滤
def double_threshold(dx_gray, dy_gray, df_gray, low, high):
	'''
	dx_gray:x方向梯度矩阵
	dy_gray:y方向梯度矩阵
	df_gray:梯度强度矩阵
	low:低阈值
	high:高阈值
	'''
    h, w = df_gray.shape
    for i in range(1, h - 1):
        for j in range(1, w - 1):
            if df_gray[i, j] < low:
                df_gray[i, j] = 0
            elif df_gray[i, j] >= high:
                df_gray[i, j] = 1
                if dy_gray[i-1, j-1] * dx_gray[i-1, j-1] > 0:  # dx，dy同向
                    if df_gray[i - 1, j + 1] > low:
                        df_gray[i - 1, j + 1] = high
                    if df_gray[i + 1, j - 1] > low:
                        df_gray[i + 1, j - 1] = high
                    if dy_gray[i-1, j-1] > dx_gray[i-1, j-1]:
                        if df_gray[i, j + 1] > low:
                            df_gray[i, j + 1] = high
                        if df_gray[i, j - 1] > low:
                            df_gray[i, j - 1] = high
                    else:
                        if df_gray[i - 1, j] > low:
                            df_gray[i - 1, j] = high
                        if df_gray[i + 1, j] > low:
                            df_gray[i + 1, j] = high
                else:
                    if df_gray[i - 1, j - 1] > low:
                        df_gray[i - 1, j - 1] = high
                    if df_gray[i + 1, j + 1] > low:
                        df_gray[i + 1, j + 1] = high
                    if math.fabs(dy_gray[i-1, j-1]) > math.fabs(dx_gray[i-1, j-1]):
                        if df_gray[i, j + 1] > low:
                            df_gray[i, j + 1] = high
                        if df_gray[i, j - 1] > low:
                            df_gray[i, j - 1] = high
                    else:
                        if df_gray[i - 1, j] > low:
                            df_gray[i - 1, j] = high
                        if df_gray[i + 1, j] > low:
                            df_gray[i + 1, j] = high
            else:
                df_gray[i, j] = 0
    return df_gray

优化结果：

优化后，运行时间增加了￣□￣｜｜，为46s，不过值得一提是，提取的轮廓质量更好了。原先想要的但是没能够输出的边，竟然输出了！ 而且花瓣的轮廓提取更加连续了。
如下上图为未优化前提取的轮廓，下图为优化后提取的轮廓。效果还是不错的。

总结：

基本实现了图像轮廓的提取，但是运行时间较长，代码还需要优化。

参考资料：

canny 算子python实现

一、简介

在图像处理中，边缘提取是一种基础的图像处理算法，常用于图像识别以及跟踪领域，为进一步分析和理解图像做准备，下面介绍两种不同的图像边缘检测方法。

二、边缘检测方法

本文介绍的边缘检测方法包括直接使用高斯滤波器检测和Canny边缘检测两种方法。

1.高斯滤波器提取边缘特征

高斯滤波是一种线性平滑滤波，可以用于消除图像中的高斯噪声。简言之，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。高斯滤波的具体操作是：用一个模板（或称卷积、掩模）扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值。下面用代码显示一下处理效果：

import cv2
img = cv2.imread("gaosi.jpg",0)
blurred = cv2.GaussianBlur(img,(11,11),0) #高斯矩阵的长与宽都是11，标准差为0
gaussImg = img - blurred
cv2.imshow("Image",gaussImg)
cv2.waitKey(0)

原图：

处理后：

2.Canny边缘特征提取

Canny边缘检测是一种非常流行的边缘检测算法，是由John F.Canny 在1986 年提出的，Canny边缘检测算法被很多人推崇为当今最优秀的边缘检测算法，关键步骤有两步，即图像去噪和设置阈值来确定边界。
图像去噪：由于边缘检测很容易受到噪声的影响，所以首先要使用高斯滤波器来去除噪声；
设置阈值：对于边界的确定需要我们来设置两个阈值，minv和maxv，当图像的灰度梯度高于maxv时，被认定为真的边界，图像的灰度小于minv时，则被舍弃，倘若灰度梯度介于两者之间，要看这个点是否与某个被认定为真的边界点相连，如果是，就被认为是真的边界点，反之则舍弃。如下图所示，A点高于阈值maxv，所以直接被认定真的边界点，B点和C点虽然都介于minv和maxv之间，但是C点与A点相连，所以也被看作是边界点，B点不与边界点相连，所以被舍弃。

下面是Canny边缘检测代码：

import cv2
img = cv2.imread("z2.jpg",0)
blurred = cv2.GaussianBlur(img,(11,11),0)
gaussImg = cv2.Canny(blurred, 10, 70)
cv2.imshow("Img",gaussImg)
cv2.waitKey(0)

原图：

处理后：

以上就是文章的所有内容，如果本文对你有帮助，可以给我点个赞，加个关注，谢谢！