Edge Detection — Canny

Edge Detection

边缘检测通过大幅度的减少待处理的数据量来简化图像分析，同时提取出图像中比较突出的信息，即有用的物体边界的结构信息[1]。边缘检测是目标识别，场景识别前期步骤。
边缘检测要满足如下三个标准：

• 低误差率。不能将非边缘处误解为边缘，如纹理噪声中会出现伪边缘。
• 边缘点应被很好地定位。已经定位的边缘必须尽可能的接近真实边缘，即检测算子标记为边缘点的一个点和真实的边缘的中心之间的距离应该最小。
• 单个边缘点响应。对于每一个真实的边缘点，检测算子应该只返回一个点。即真实边缘周围的局部最大数应该是最小的。尽可能避免检测出两条边距离特别小的双沟边的现象。

另外，边缘检测和定位之间存在一个测不准的原则（uncertainly principle）

• 可以通过改变检测算子的空间宽度来获得边缘检测和定位之间的平衡。

Canny 算子

Canny算子假设

图像除了有光滑变换的区域和边缘之外，还有可叠加的高斯白噪声。

• 不存在角，即canny算子不能用来检测角。
• 如果需要检测角，可以用Harris算子

Canny算子

1 通过高斯低通滤波来对输入图像进行滤波。

• 通过滤波去除图像中的噪声，以免噪声在检测时产生伪边缘。

2 计算梯度幅度图像和角度图像

• 令𝑓(𝑥,𝑦)为经高斯滤波后得到的平滑图像
• 则其梯度幅度图像为：
• 角度图像为：
  
  

3 对梯度幅度图像应用非极大抑制 ( non-maximum suppression )。

• 沿着检测出的边缘的方向，对比每一个像素点的梯度幅度和沿着梯度方向（边缘法线）前后两个点梯度幅度的大小，如果两点梯度幅度值都比中间点小，则该点可能为边缘点将其保留。否则该点不是边缘点将其抑制，后续不再考虑。

4 使用双阈值处理和连通性分析来检测与连接边缘。

• 首先设置一高一低两个阈值，跟踪未被抑制的像素。如果像素的梯度幅度低于较低的阈值，则将像素置零。如果像素的梯度幅度高于较高的阈值，则保留该像素点。如果像素梯度幅度介于两个阈值之间，则若与该像素8连通的点有边缘点，则保留该像素点，否则置零（该点可能为噪点）。

使用Canny算子进行边缘检测

Step 1: Gaussian Filter

下图是一个滤波器：

Step 2 ：Sobel Operator

• 通过Sobel算子来计算图像中每一个像素点的二维空间梯度(2-D spatial gradient ) 。
• 在数字图像处理中，使用差分来代替高数中的求导。
  在每个像素点得到其梯度分量𝐺_𝑥和𝐺_y后，通过使用绝对值来近似梯度幅度：
  • |𝐺|=|𝐺_𝑥 |+|𝐺_𝑦 |
• 边的方向，即角度为：
  • 𝜃=arctan⁡(𝐺_y∕𝐺_x )
    当𝐺_x=0且𝐺_y不为0时，𝜃为90度。

下图是一个Sobel算子：

Step 3 ：边缘检测

• 梯度幅度图像通常在局部极大值附近包含一些宽脊，为了细化宽脊要使用非极大值抑制

• 水平方向

• 45°方向 （沿着正对角线）

• 垂直方向

• -45°方向（沿着负对角线）

• 由边缘法线的方向确定边缘的方向

  • 边缘法线方向即𝜃

下图代表法线对应相应的边缘范围：

Step 4 ：非极大值抑制

令𝑑_1 、 𝑑_2 、 𝑑_3 、𝑑_4分别代表上述四个基本方向：水平方向，45°，垂直方向，-45°，非极大值抑制方案如下：

• 寻找最接近𝜃(𝑥,𝑦)的方向𝑑_𝑘 。
• 令K表示‖∇𝑓(𝑥, 𝑦)‖在(𝑥,𝑦)处的值。若K或大于𝑑_𝑘梯度方向上点(𝑥,𝑦)的前后两个邻点处的‖∇𝑓(𝑥, 𝑦)‖值，则令该点像素值为K，否则为零（抑制）。

Step 4 ：双阈值处理

• 使用一个低阈值𝑇_𝐿和高阈值𝑇_𝐻。实验^([1])表明高低阈值比率应在2:1到3:1的范围内。
• 如果像素点(𝑥, 𝑦)的梯度幅度‖∇𝑓(𝑥, 𝑦)‖ ≤ 𝑇_𝐿，则该像素点不是边缘点，并令其像素值为0
• 如果像素点(𝑥, 𝑦)的梯度幅度‖∇𝑓(𝑥, 𝑦)‖ ≥𝑇_𝐻，则该像素点是边缘点。
• 如果像素点(𝑥, 𝑦)的梯度幅度‖∇𝑓(𝑥, 𝑦)‖的值在高低阈值之间，则看该像素点的8邻域中的点有没有边缘点，如果有则该像素点是边缘点；否则不是，将其像素值置0 。

代码&结果展示

相应代码如下：
Python3.8 实现

import cv2 as cv
import numpy as np

class Canny:
    # 将彩色图像转化为灰度图像
    def BGR2GRAY(self, src):
        b = src[:, :, 0].copy()
        g = src[:, :, 1].copy()
        r = src[:, :, 2].copy()

        out = 0.2126*r + 0.7152*g + 0.0722*g
        out = out.astype(np.uint8)

        return out

    # 对灰度图像进行高斯滤波
    def gaussion_filter(self, img, k_size=3, sigma=1.3):
        if len(img.shape) == 3:
            H, W, C = img.shape
            gray = False
        else:
            img = np.expand_dims(img, axis=-1)
            H, W, C = img.shape
            gray = True

            ## Zero padding
        pad = k_size // 2
        out = np.zeros([H + pad * 2, W + pad * 2, C], dtype=np.float)
        out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float)

        ## prepare Kernel
        K = np.zeros((k_size, k_size), dtype=np.float)
        for x in range(-pad, -pad + k_size):
            for y in range(-pad, -pad + k_size):
                K[y + pad, x + pad] = np.exp(- (x ** 2 + y ** 2) / (2 * sigma * sigma))
        # K /= (sigma * np.sqrt(2 * np.pi))
        K /= (2 * np.pi * sigma * sigma)
        K /= K.sum()

        tmp = out.copy()

        # 滤波
        for y in range(H):
            for x in range(W):
                for c in range(C):
                    out[pad + y, pad + x, c] = np.sum(K * tmp[y: y + k_size, x: x + k_size, c])

        out = np.clip(out, 0, 255)
        out = out[pad: pad + H, pad: pad + W]
        out = out.astype(np.uint8)

        if gray:
            out = out[..., 0]

        return out
    # sobel 滤波
    def sobel_filter(self, src, k_size=3):
        if len(src.shape) == 3:
            H, W, C = src.shape
        else:
            H, W = src.shape

        ## 零填充
        pad = k_size // 2
        out = np.zeros((H+pad*2, W+pad*2), dtype=np.float)
        out[pad:pad+H, pad:pad+W] = src.copy().astype(np.float)
        tmp = out.copy()

        out_v = out.copy()
        out_h = out.copy()

        ## Sobel vertical
        Kv = [[1., 2., 1.], [0., 0., 0.], [-1., -2., -1.]]
        ## Sobel horizontal
        Kh = [[1., 0., -1.], [2., 0., -2.], [1., 0., -1.]]

        # 滤波
        for y in range(H):
            for x in range(W):
                out_v[pad + y, pad + x] = np.sum(Kv * (tmp[y: y + k_size, x: x + k_size]))
                out_h[pad + y, pad + x] = np.sum(Kh * (tmp[y: y + k_size, x: x + k_size]))

        out_v = np.clip(out_v, 0, 255)
        out_h = np.clip(out_h, 0, 255)

        out_v = out_v[pad: pad + H, pad: pad + W]
        out_v = out_v.astype(np.uint8)
        out_h = out_h[pad: pad + H, pad: pad + W]
        out_h = out_h.astype(np.uint8)

        return out_v, out_h

    def get_edge_angle(self, fx, fy):
        # get edge strength
        edge = np.sqrt(np.power(fx.astype(np.float32), 2) + np.power(fy.astype(np.float32), 2))
        edge = np.clip(edge, 0, 255) #将像素值固定在0-255之间

        fx = np.maximum(fx, 1e-10) # 避免fx=0
        # fx[np.abs(fx) <= 1e-5] = 1e-5

        # get edge angle
        angle = np.arctan(fy / fx)

        return edge, angle

    # 根据角度范围确定边缘方向
    def angle_quantization(self, angle):
        angle = angle / np.pi * 180
        angle[angle < -22.5] = 180 + angle[angle < -22.5]
        _angle = np.zeros_like(angle, dtype=np.uint8)
        _angle[np.where(angle <= 22.5)] = 0
        _angle[np.where((angle > 22.5) & (angle <= 67.5))] = 45
        _angle[np.where((angle > 67.5) & (angle <= 112.5))] = 90
        _angle[np.where((angle > 112.5) & (angle <= 157.5))] = 135

        return _angle

    # 非极大值抑制
    def non_maximum_suppression(self, angle, edge):
        H, W = angle.shape
        _edge = edge.copy()

        for y in range(H):
            for x in range(W):
                if angle[y, x] == 0:
                    dx1, dy1, dx2, dy2 = -1, 0, 1, 0
                elif angle[y, x] == 45:
                    dx1, dy1, dx2, dy2 = -1, 1, 1, -1
                elif angle[y, x] == 90:
                    dx1, dy1, dx2, dy2 = 0, -1, 0, 1
                elif angle[y, x] == 135:
                    dx1, dy1, dx2, dy2 = -1, -1, 1, 1
                if x == 0:
                    dx1 = max(dx1, 0)
                    dx2 = max(dx2, 0)
                if x == W - 1:
                    dx1 = min(dx1, 0)
                    dx2 = min(dx2, 0)
                if y == 0:
                    dy1 = max(dy1, 0)
                    dy2 = max(dy2, 0)
                if y == H - 1:
                    dy1 = min(dy1, 0)
                    dy2 = min(dy2, 0)
                if max(max(edge[y, x], edge[y + dy1, x + dx1]), edge[y + dy2, x + dx2]) != edge[y, x]:
                    _edge[y, x] = 0

        return _edge

    # 使用双阈值处理和连通性分析来检测边缘
    def hysterisis(self, edge, HT=100, LT=30):
        H, W = edge.shape

        # 双阈值
        edge[edge >= HT] = 255
        edge[edge <= LT] = 0

        _edge = np.zeros((H + 2, W + 2), dtype=np.float32)
        _edge[1: H + 1, 1: W + 1] = edge

        ## 8 - 邻域
        nn = np.array(((1., 1., 1.), (1., 0., 1.), (1., 1., 1.)), dtype=np.float32)

        # 如果梯度幅度值在LT和LH之间
        for y in range(1, H + 2):
            for x in range(1, W + 2):
                if _edge[y, x] < LT or _edge[y, x] > HT:
                    continue
                if np.max(_edge[y - 1:y + 2, x - 1:x + 2] * nn) >= HT: #判断八连通区域中是否有边缘点
                    _edge[y, x] = 255
                else:
                    _edge[y, x] = 0

        edge = _edge[1:H + 1, 1:W + 1]

        return edge

def main():
    src = cv.imread("C:/Users/Odysseus96/Pictures/Image/001.JPG").astype(np.float32)
    canny = Canny()
    # 转变成灰度图像
    out = canny.BGR2GRAY(src)
    cv.imshow("Input", out)
    # 高斯滤波
    gaussian = canny.gaussion_filter(out)
    # sobel 滤波，获得gx， gy
    fx, fy = canny.sobel_filter(gaussian, 3)
    # 计算梯度幅度图像和角度图像
    edge, angle = canny.get_edge_angle(fx, fy)

    # cv.imshow("output", edge)
    # cv.imwrite("C:/Users/Odysseus96/Pictures/Image/edge.JPG", edge)
    # cv.waitKey(0)
    # cv.destroyAllWindows()

    # 根据角度范围确定边缘方向
    angle = canny.angle_quantization(angle)
    # 对梯度幅度图像应用非极大抑制
    edge = canny.non_maximum_suppression(angle, edge)
    # 使用双阈值处理和连通性分析来检测与连接边缘
    out = canny.hysterisis(edge, 50, 20).astype(np.uint8)

    # cv.imwrite("C:/Users/Odysseus96/Pictures/Image/canny_001.JPG", out)

    cv.imshow("output", out)
    cv.waitKey(0)
    cv.destroyAllWindows()

if __name__ == '__main__':
    main()

对于灰度图像：

运行结果如下：
检测到的边缘图像

参考文献：

[1] J. Canny, “A Computational Approach to Edge Detection,” IEEE Trans. Pattern Analysis and Machine Intelligence, 8(6), pp. 679-698, 1986.
[2] digital image processing 4th edition