canny算法

import cv2 as cv
import numpy as np
import math
filepath = 'C:\\Users\\10107472\\Desktop\\tradition_cv\\image\\0.jpg'
# step one 灰度化
"""
计算公式：
Gray(i,j) = [R(i,j) + G(i,j) + B(i,j)] / 3
or :
Gray(i,j) = 0.299 * R(i,j) + 0.587 * G(i,j) + 0.114 * B(i,j)
"""
# 读取图片
img = cv.imread('0.jpg')  # 读取文件位置图片，此处修改图片位置
# BGR 转换成 RGB 格式
img_rgb = cv.cvtColor(img, cv.COLOR_BGR2RGB)
# 灰度化
img_gray = np.dot(img_rgb[..., :3], [0.299, 0.587, 0.114])
print('灰度后图像的形状 = ', img_gray.shape)                                  # 查看灰度图像形状


# step two 去除噪音（ 使用 5x5 的高斯滤波器 ）
# 生成高斯滤波器
"""
要生成一个 (2k+1)x(2k+1) 的高斯滤波器，滤波器的各个元素计算公式如下：
H[i, j] = (1/(2*pi*sigma**2))*exp(-1/2*sigma**2((i-k-1)**2 + (j-k-1)**2))
"""
sigma1 = sigma2 = 1.4
gau_sum = 0
gaussian = np.zeros([5, 5])
for i in range(5):
    for j in range(5):
        gaussian[i, j] = math.exp((-1/(2*sigma1*sigma2))*(np.square(i-2-1) + np.square(j-2-1)))/(2*math.pi*sigma1*sigma2)
        gau_sum = gau_sum + gaussian[i, j]
gaussian = gaussian / gau_sum  # 归一化处理
# 高斯滤波
W, H = img_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(img_gray[i:i+5, j:j+5] * gaussian)

print('高斯后的图像形状 = ', new_gray.shape)                                  # 查看高斯滤波后图像的形状

# step three 计算梯度幅值

"""
:type: image which after smooth
:rtype:
    dx: gradient in the x direction
    dy: gradient in the y direction
    M: gradient magnitude
    theta: gradient direction
"""

W, H = new_gray.shape
dx = np.zeros([W-1, H-1])
dy = np.zeros([W-1, H-1])
M = np.zeros([W-1, H-1])
theta = np.zeros([W-1, H-1])
for i in range(W-1):
    for j in range(H-1):
        dx[i, j] = new_gray[i+1, j] - new_gray[i, j]   # 计算x的梯度
        dy[i, j] = new_gray[i, j+1] - new_gray[i, j]
        M[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j]))    # 计算图像梯度幅值作为图像强度值
        theta[i, j] = math.atan(dx[i, j] / (dy[i, j] + 0.000000001))  # 计算  θ - artan(dx/dy) ，计算梯度方向
print('x,y方向上的梯度,幅值与角度 = ', dx.shape)
cv.imshow('dx ', dx)
cv.imshow(' dy ', dy)
cv.imshow('M ', M)
cv.imshow('theta', theta)
# step four NMS
d = np.copy(M)    # M为幅度值
W, H = M.shape
NMS = np.copy(d)
NMS[0, :] = NMS[W-1, :] = NMS[:, 0] = NMS[:, H-1] = 0    # 表示四周取0，即第一行、最后一行、第一列、最后一列
for i in range(1, W-1):
    for j in range(1, H-1):
        # 如果当前梯度为0，该点就不是边缘点，原因在于边缘点像素与周围像素变化很大，因此导致梯度至少不为0
        if M[i, j] == 0:
            NMS[i, j] = 0
        else:
            gradX = dx[i, j]        # 当前点 x 方向导数
            gradY = dy[i, j]        # 当前点 y 方向导数
            gradTemp = d[i, j]      # 当前梯度点
            if np.abs(gradY) > np.abs(gradX):               # 如果 y 方向梯度值比较大，说明导数方向趋向于 y 分量
                weight = np.abs(gradX) / np.abs(gradY)      # 权重
                grad2 = d[i-1, j]
                grad4 = d[i+1, j]
                if gradX * gradY > 0:    # 如果 x, y 方向导数符号一致
                    grad1 = d[i-1, j-1]
                    grad3 = d[i+1, j+1]
                else:      # 如果 x，y 方向导数符号相反
                    grad1 = d[i-1, j+1]
                    grad3 = d[i+1, j-1]
            else:       # 如果 x 方向梯度值比较大
                weight = np.abs(gradY) / np.abs(gradX)
                grad2 = d[i, j-1]
                grad4 = d[i, j+1]          
                if gradX * gradY > 0:       # 如果 x, y 方向导数符号一致
                    grad1 = d[i+1, j-1]
                    grad3 = d[i-1, j+1]
                else:        # 如果 x，y 方向导数符号相反
                    grad1 = d[i-1, j-1]
                    grad3 = d[i+1, j+1]
            # 利用 grad1-grad4 对梯度进行插值
            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               # 不可能是边缘点
cv.imshow('image after nms', NMS)
W, H = NMS.shape
DT = np.zeros([W, H])
TL = 0.1 * np.max(NMS)      # 定义低阈值
TH = 0.3 * np.max(NMS)      # 定义高阈值
for i in range(1, W-1):
    for j in range(1, H-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
cv.imshow('final image', DT)
cv.waitKey( )

NMS 是要找出局部最大值，因此，需要将当前的像素的梯度，与其他方向进行比较。如下图所示，g1，g2，g3，g4 分别是 C 八个领域中的 4 个点，蓝线是 C 的梯度方向。如果 C 是局部最大值的话，C 点的梯度幅值就要大于梯度方向直线与 g1g2，g4g3 两个交点的梯度幅值，即大于点 dTemp1 和 dTemp2 的梯度幅值。上面提到这种方法无法达到最好的效果，因为 dTemp1 和 dTemp2 不是整像素，而是亚像素。亚像素的意思就是在两个物理像素之间还有像素。那么，亚像素的梯度幅值怎么求？可以使用线性插值的方法，计算 dTemp1 在 g1，g2 之间的权重，就可以得到其梯度幅值。计算公式如下：

weight = |gx| / |gy| or |gy| / |gx|
dTemp1 = weight*g1 + (1-weight)*g2
dTemp2 = weight*g3 + (1-weight)*g4

  计算时分两种情况（都是比较当前像素与dtemp1与dtemp2的大小，大于这两个值则保留，小于其中任意一个则将其值取0）：1）下面两幅图是 y 方向梯度值比较大的情况，即梯度方向靠近 y 轴。所以，g2 和 g4 在 C 的上下位置，此时 weight = |gy| / |gx| 。左边的图是 x，y 方向梯度符号相同的情况，右边是 x，y 方向梯度符号相反的情况。

下面两幅图是 x 方向梯度值比较大的情况，即梯度方向靠近 x 轴。所以，g2 和 g4 在 C 的左右位置，此时 weight = |gy| / |gx| 。左边的图是 x，y 方向梯度符号相同的情况，右边是 x，y 方向梯度符号相反的情况。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

最终结果：