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.效果
校医院还是要检测一个的。