赞
踩
Canny算子是一种较为好用的边缘检测算子,边缘检测效果良好。利用Python实现,原速度约为5分钟以上,不能够满足使用需要。
本文利用Python实现,并注意使用numba支持的数据类型及写法进行,故而可以利用numba进行加速,对5000*3000像素图像提取边缘,整体流程用时13秒,满足使用需要。
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
用于提供一个高斯滤波核
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
暂时只实现了两种梯度,后面可能会增加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
# -*- 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
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则不是。
校医院还是要检测一个的。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。