Canny算子的Python实现

1.引言

Canny算子是一种较为好用的边缘检测算子，边缘检测效果良好。利用Python实现，原速度约为5分钟以上，不能够满足使用需要。

本文利用Python实现，并注意使用numba支持的数据类型及写法进行，故而可以利用numba进行加速，对5000*3000像素图像提取边缘，整体流程用时13秒，满足使用需要。

2.预备功能函数

1.卷积函数

def convolution(image, kernel):
    """
    This function is used for convolution.
    Args:
        image: image input.
        kerbel: the filter.
    """

    kernel_height, kernel_width = kernel.shape
    height, width = image.shape
    kernel_size = kernel_height

    half_size = np.floor(kernel_size/2)

    

    @numba.jit(nopython = True)
    def conv(image, kernel, half_size, height, width):

        result = np.zeros((height, width), dtype = np.float32)
        for row in range(half_size, height - half_size):
            for col in range(half_size, width - half_size):
                sum_var = 0
                for v_row in range(-1*half_size, half_size+1):
                    for v_col in range(-1*half_size, half_size+1):
                        sum_var = sum_var + image[row+v_row, col+v_col] * kernel[v_row, v_col]

                result[row, col] = sum_var

        return result

    result = conv(image, kernel, half_size, height, width)
    return result

2.高斯模板函数

用于提供一个高斯滤波核

def generate_Guassian_template(kernel_size = 3, sigma = 1):
                
    template = np.zeros((kernel_size, kernel_size), \
                    dtype = np.float32)
    halfsize = np.floor(kernel_size/2).astype(np.int16)

    @numba.jit(nopython = True)
    def gaussian2d(x, y, sigma):
        
        result = np.exp(-(np.power(x, 2)+np.power(y, 2))/(2*np.power(sigma, 2)))
        return result

    @numba.jit(nopython = True)
    def generate(template, halfsize, sigma):

        for v_row in range(-1*halfsize, halfsize+1):
            for v_col in range(-1*halfsize, halfsize+1):
                template[v_row+halfsize, v_col+halfsize] = gaussian2d(v_row, v_col, sigma)

        element_sum = np.sum(template)
        template = template/element_sum
        
        return template

    re = generate(template, halfsize, sigma)
    return re

3. 梯度类

暂时只实现了两种梯度，后面可能会增加Sobel梯度。

# -*- coding: utf-8 -*-
# @author: Dasheng Fan
# @email: fandasheng1999@163.com

import numpy as np
import numba
import sys

sys.path.append('../')
from utils.gray_processing import gray_processing

class RobertsGradient(object):
    """
    Roberts is a kind of gradient of the image.

    """

    def __init__(self):
        """
        """
        
    def calculate(self, image):
        """
        Culculate the Roberts gradient gx and gy.

        Args:
            image: image to input.

        Return:
            gx, gy
        """

        if len(image.shape) != 2:
            image = gray_processing(image)

        @numba.jit(nopython = True)
        def run(image):

            height, width = image.shape

            gradient_gy_img = np.zeros((height, width), dtype = np.float32)
            gradient_gx_img = np.zeros((height, width), dtype = np.float32)

            for row in range(1, height-1):
                for col in range(1, width-1):
                    gradient_gx_img[row, col] = image[row, col] - image[row+1, col+1]
                    gradient_gy_img[row, col] = image[row+1, col] - image[row, col+1]

            return gradient_gx_img, gradient_gy_img

        gx, gy = run(image)

        return gx, gy


class Gradient(object):

    def __init__(self):
        pass
    
    
    def calculate(self, image):

        @numba.jit(nopython=True)
        def run(image):
            height, width = image.shape

            gx = np.zeros((height, width), dtype = np.float32)
            gy = np.zeros((height, width), dtype = np.float32)

            for row in range(1, height-1):
                for col in range(1, width-1):
                    gx[row, col] = image[row, col+1] - image[row, col]
                    gy[row, col] = image[row+1, col] - image[row, col]

            return gx, gy

        gx, gy = run(image)
        return gx, gy

4.灰度化函数

# -*- coding: utf-8 -*-
# @author: Dasheng Fan
# @email: fandasheng1999@163.com

"""
Gray processing.
"""

import numpy as np
from matplotlib import pyplot as plt
import sys
import os

def gray_processing(image):
    """
    For gray processing.
    
    Args:
        image: A 3 channels image.
    
    Reture:
        A gray image.

    Reference:
    "https://baike.baidu.com/item/%E7%81%B0%E5%BA%A6%E5%8C%96/3206969?fr=aladdin"
    """

    if len(image.shape) != 3:
        raise Exception("The channel is wrong.")

    u = np.power(image[:, :, 0], 2.2) + np.power(1.5*image[:, :, 1], 2.2) + \
        np.power(0.6*image[:, :, 2], 2.2)

    d = 1 + np.power(1.5, 2.2) + np.power(0.6, 2.2)

    gray = np.power(u/d, 1/2.2)

    return gray

if __name__ == '__main__':
    pass

3.实现

1.高斯滤波，参考：
https://www.jianshu.com/p/73e6ccbd8f3f
2.梯度计算：无参考
3.极大值抑制（Canny的很复杂）
主要参考：https://blog.csdn.net/kezunhai/article/details/11620357
4.双阈值检测，抑制孤立低阈值点
参考：https://www.cnblogs.com/techyan1990/p/7291771.html

代码：

# -*- coding: utf-8 -*-
# @author: Dasheng Fan
# @email: fandasheng1999@163.com

"""
Canny operator.
"""
import numpy as np
import numba
from operators.diptools import generate_Guassian_template, convolution
from operators.gradient import RobertsGradient,Gradient
from utils.gray_processing import gray_processing
from utils.display import display, doubledisplay


class Canny(object):
    """
    Canny operator is used for line detecting.
    """

    def __init__(self, 
            Guassian_kernel_size = 3, 
            Guassian_sigma = 1, 
            high_threshhold_ratio = 0.15,
            low_threshhold_ratio = 0.08):

        self._G_kernel_size = Guassian_kernel_size
        self._G_sigma = Guassian_sigma
        self._h_thresh_r = high_threshhold_ratio
        self._l_thresh_r = low_threshhold_ratio

    def detect(self, image):

        if len(image.shape) != 2:
            image = gray_processing(image)

        height, width = image.shape

        g_template = generate_Guassian_template()
        image = convolution(image, g_template)

        gradient = Gradient()
        gx, gy = gradient.calculate(image)

        @numba.jit(nopython = True)
        def nonmaxsuppress(gx,gy):
            """
            Reference:
            https://blog.csdn.net/kezunhai/article/details/11620357
            """
            
            g = np.sqrt(np.add(np.power(gx, 2), np.power(gy, 2)))
            angle = np.arctan(gy/gx)
            
            height, width = g.shape
            flags = np.zeros((height, width), dtype = np.float32)

            for row in range(1, height-1):
                for col in range(1, width - 1):

                    local_g = g[row, col]

                    if np.abs(gy[row, col])>np.abs(gx[row, col]):
                        
                        if gy[row, col] == 0:
                            weight = 1
                        else:
                            weight = np.abs(gx[row,col]/gy[row, col])
                        g2 = g[row-1, col]
                        g4 = g[row+1, col]
                        
                        if np.multiply(gx[row,col], gy[row, col])>0:
                            g1 = g[row-1, col-1]
                            g3 = g[row+1, col+1]
                        else:
                            g1 = g[row-1, col+1]
                            g3 = g[row+1, col-1]
                    else:
                        if gx[row, col] == 0:
                            weight = 1
                        else:
                            weight = np.abs(gy[row, col]/gx[row, col])
                        g2 = g[row, col-1]
                        g4 = g[row, col+1]

                        if np.multiply(gx[row, col],gy[row,col])>0:
                            g1 = g[row-1, col-1]
                            g3 = g[row+1, col+1]
                        else:
                            g1 = g[row+1, col-1]
                            g3 = g[row-1, col+1]

                    inter_g1 = weight*g1 + (1-weight)*g2
                    inter_g2 = weight*g3 + (1-weight)*g4 
                    local_g = g[row, col]

                    if local_g >=inter_g1 and local_g>=inter_g2:
                        flags[row, col] = 1
            
            return flags

        flags1 = nonmaxsuppress(gx, gy)
        
        @numba.jit(nopython = True)
        def _double_threshhold_suppress(
                                high_threshhold_ratio, 
                                low_threshhold_ratio, 
                                gx, gy):
            """
            The theory can be found here:
            https://www.cnblogs.com/techyan1990/p/7291771.html
            Only the theory is refered, the codes are different.
            """

            g = np.sqrt(np.add(np.power(gx, 2), np.power(gy, 2)))
            height, width = g.shape
            max_g = np.max(g)
            high_thresh = max_g * high_threshhold_ratio
            low_thresh = max_g * low_threshhold_ratio

            flags = np.zeros((height, width), dtype = np.float32)
            
            for row in range(1, height-1):
                for col in range(1, width - 1):
                    
                    if g[row, col] >= high_thresh:
                        flags[row, col] = 1
                    elif g[row, col] > low_thresh and g[row, col] < high_thresh:
                        # 不洋嚯了，写汉语
                        # 这里检查一下八邻域， 如果有强边缘就认为弱边缘是边缘点
                        for var_y in range(-1, 2):
                            for var_x in range(-1, 2):
                                if g[row+var_y, row+var_x] > high_thresh:
                                    flags[row, col] = 1
                                    break
                    else:
                        flags[row, col] = 0
            
            return flags

        flags2 = _double_threshhold_suppress(self._h_thresh_r, self._l_thresh_r, gx, gy)

        flags = np.multiply(flags1, flags2)

        return flags

返回的flags是标志矩阵，矩阵中为1的元素对应的点是边缘点，为0则不是。

4.效果

校医院还是要检测一个的。