canny算法原理及实现

一、canny算法原理

1. 灰度化
2. 高斯滤波

•  通过对待滤波像素及其相邻点进行加权均值计算，可以去掉图像上的噪点 

1. 提取边缘
2. 非极大值抑制

• 梯度值较大的点可能是真正的边缘像素点，也可能是由颜色变化引起或噪点，想要过滤掉非边缘可以使用非极大值抑制法
• 像素点和沿梯度方向上相邻的两个像素点进行比较，如果此像素点是三个点中梯度值最大的，则保留，否则置为0

1. 双阈值法

• 给定两个阈值，如果像素点的梯度值高于高阈值，则为强边缘，如果低于低阈值，则不是边缘，介于低阈值和高阈值之间的为弱边缘，对于弱边缘又分两种情况：一种是真实边缘附近的点，一种孤立的点，不是真实边缘
• 对于弱边缘，需要判断此点与其相邻的八个点中是否有强边缘，如果有，保留此点，如果没有，则删除

二、canny算法实现

#coding=utf-8

import numpy as np
import cv2


class Canny():
    def img2gray(self, img):
        b = img[:, :, 0].copy()
        g = img[:, :, 1].copy()
        r = img[:, :, 2].copy()
        out = b*0.0722 + g*0.7152 + r*0.2126
        out = out.astype(np.uint8)
        return out

    def gussian_filter(self, img, ksize=3, sigma=1.4):
        H, W= img.shape
        pad = ksize // 2
        out = np.zeros([H+2*pad, W+2*pad], dtype=np.float)
        out[pad:pad+H, pad:pad+W] = img.copy().astype(np.float)

        K = np.zeros((ksize, ksize), dtype=np.float)

        for x in range(-pad, -pad+ksize):
            for y in range(-pad, -pad+ksize):
                K[y+pad, x+pad] = np.exp(-(x**2 + y**2)/2*sigma*sigma)

        K /= (2 * np.pi*sigma*sigma)
        K /= K.sum()

        tmp = out.copy()
        for y in range(H):
            for x in range(W):
                out[y + pad, x + pad] = np.sum(K*tmp[y: y+ksize, x :x+ksize])

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

    def sobel_filter(self, img, ksize=3):
        pad = ksize // 2
        H, W= img.shape
        out = np.zeros([H + pad * 2, W + pad * 2], dtype=np.float)
        out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float)
        tmp = out.copy()

        out_y = out.copy()
        out_x = out.copy()

        kx = [[1., 0., -1.], [2., 0., -2.], [1., 0., -1.]]
        ky = [[1., 2., 1.], [0., 0., 0.], [-1., -2., -1.]]

        for y in range(H):
            for x in range(W):
                out_y[pad+y, pad+x] = np.sum(ky*tmp[y:y+ksize, x:x+ksize])
                out_x[pad+y, pad+x] = np.sum(kx*tmp[y:y+ksize, x:x+ksize])

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

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

        return out_x, out_y

    def get_angle(self, out_x, out_y):
        edge = np.sqrt(np.power(out_x.astype(np.float32), 2) + np.power(out_y.astype(np.float32), 2))
        edge = np.clip(edge, 0, 255)
        out_x = np.maximum(out_x, 1e-10)
        angle = np.arctan(out_y/out_x)
        return edge, angle

    def angle_handle(self, angle):
        angle = angle / np.pi * 180
        angle[angle < -22.5] = 180 + angle[angle < -22.5]
        new_angle = np.zeros_like(angle, dtype=np.uint8)

        new_angle[np.where(angle <= 22.5)] = 0
        new_angle[np.where((angle > 22.5) & (angle <= 67.5))] = 45
        new_angle[np.where((angle > 67.5) & (angle <= 112.5))] = 90
        new_angle[np.where((angle > 112.5) & (angle <= 157.5))] = 135

        return new_angle

    def non_max_sus(self, edge, angle):
        H, W = edge.shape
        new_edge = edge.copy()
        for y in range(0, H):
            for x in range(0, W):
                if angle[y, x] == 0:
                    x1, y1, x2, y2 = -1, 0, 1, 0
                elif angle[y, x] == 45:
                    x1, y1, x2, y2 = -1, -1, 1, 1
                elif angle[y, x] == 90:
                    x1, y1, x2, y2 = 0, -1, 0, 1
                elif angle[y, x] == 135:
                    x1, y1, x2, y2 = -1, 1, 1, -1,

                if x == 0:
                    x1 = max(x1, 0)
                    x2 = max(x2, 0)
                if y == 0:
                    y1 = max(y1, 0)
                    y2 = max(y2, 0)
                if x == W-1:
                    x1 = min(x1, 0)
                    x2 = min(x2, 0)
                if y == H - 1:
                    y1 = min(y1, 0)
                    y2 = min(y2, 0)

                if max(max(edge[y, x], edge[y + y1, x + x1]), edge[y + y2, x + x2]) != edge[y, x]:
                    new_edge[y, x] = 0

        return new_edge

    def hysterisis(self, edge, HT = 100, LT = 10):
        H, W = edge.shape

        edge[np.where(edge >= HT)] = 255
        edge[np.where(edge <= LT)] = 0

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

        K = np.array(((1, 1, 1),(1, 0, 1),(1, 1, 1)))

        for y in range(0, H + 2):
            for x in range(0, W + 2):
                if new_edge[y, x] < LT or new_edge[y, x ]> HT:
                    continue
                if np.max(new_edge[y-1: y+2, x - 1: x + 2] * K) >= HT:
                    new_edge[y, x] = 255
                else:
                    new_edge[y, x] = 0
        edge = new_edge[1: H:1, 1: W+1]

        return edge

    def execute_func(self, img):
        gray = self.img2gray(img)
        gaussian = self.gussian_filter(gray)
        x, y = self.sobel_filter(gaussian)
        edge, angle = self.get_angle(x, y)
        new_angle = self.angle_handle(angle)
        new_edge = self.non_max_sus(edge, new_angle)
        edge_out = self.hysterisis(new_edge)
        return edge_out

test = Canny()
img = cv2.imread('lenna.png')
edge = test.execute_func(img)
edge = edge.astype(np.uint8)

cv2.imshow('img', edge)
cv2.waitKey(0)
cv2.destroyAllWindows()