门户网站的设计,深圳网站建设公司服务怎么做,网站编程图,郑州大型网站建设价格Canny边缘检测原理及其python实现 转载请注明出处#xff1a;️ Sylvan Ding Canny边缘检测算法
多数分割算法均基于灰度值的两个基本性质之一#xff1a;不连续性和相似性。图像的边缘检测就是基于灰度的不连续性#xff08;灰度突变#xff09;来实现图像分割的。…Canny边缘检测原理及其python实现 转载请注明出处©️ Sylvan Ding Canny边缘检测算法
多数分割算法均基于灰度值的两个基本性质之一不连续性和相似性。图像的边缘检测就是基于灰度的不连续性灰度突变来实现图像分割的。
边缘检测器是设计用来检测边缘像素的局部图像处理方法。一条线可视为一条边缘线段该线两侧的背景灰度要么远亮于该线像素的灰度要么远暗于该线像素的灰度这样的线称之为“屋顶边缘”。
虽然Canny[1968]算法复杂但其边缘检测的表现却相当优秀。
工具介绍
为了方便计算我们将所有灰度都标定到 [0,1]归一化这样就无需考虑系数的问题了。
二阶偏微分
∂f∂xf(x1)−f(x)\frac{\partial f}{\partial x}f(x1)-f(x)∂x∂ff(x1)−f(x)
∂2f∂x2f(x1)f(x−1)−2f(x)\frac{\partial^{2} f}{\partial x^{2}}f(x1)f(x-1)-2 f(x)∂x2∂2ff(x1)f(x−1)−2f(x)
在实验3: 拉普拉斯锐化中曾讨论过用差分定义的导数来表征和检测灰度的突变同时得到了如下结论
二阶导数仅在“灰度斜坡”的开始、结尾处不为零一阶导数在整个“灰度斜坡”上均不为零因此一阶导数产生的边缘粗于二阶导数同时二阶导数变化的激烈度明显高于一阶导数因此二阶导数更能增强细节包括噪声、孤立点、边缘等在“斜坡”和“台阶”边缘二阶导数的符号相反导致了“双边缘效应”也可用于确定边缘是从亮到暗过渡还是从暗到亮过渡
对于一条边缘来说二阶导数具有如下性质
对图像中的边缘二阶导数生成两个值双边缘效应二阶导数的零交叉点可用于定位粗边缘中心 高斯低通滤波
不难证明微弱的可见噪声会对边缘检测的上述两个关键导数产生严重影响特别是对二阶偏导数的影响非常剧烈。所以在执行边缘检测前需要对图像进行平滑处理降噪。
高斯滤波是一种线性平滑滤波适用于消除高斯噪声广泛应用于图像处理的减噪过程。通俗的讲高斯滤波就是对整幅图像进行加权平均的过程每一个像素点的值都由其本身和邻域内的其他像素值经过加权平均后得到。
二维高斯函数μ0\mu0μ0, 中心点为原点的定义如下
G(x,y)e−x2y22σ2G(x, y)\mathrm{e}^{-\frac{x^{2}y^{2}}{2 \sigma^{2}}}G(x,y)e−2σ2x2y2 注G(x,y)G(x,y)G(x,y) 的系数 12πσ2\frac{1}{2\pi\sigma^2}2πσ21 可忽略因为在计算得到权重矩阵后还要对这个矩阵进行归一化。 用 GGG 和输入图像 fff 卷积形成平滑后的图像 fsf_sfs
fs(x,y)G(x,y)∗f(x,y)f_s(x,y)G(x,y)\ast f(x,y)fs(x,y)G(x,y)∗f(x,y)
梯度和梯度方向
为了在 fff 中寻找边缘的强度和方向所以引入梯度 ∇f\nabla f∇f定义如下
∇f≡grad(f)≡[gxgy][∂f∂x∂f∂y]\nabla f \equiv \operatorname{grad}(f) \equiv\left[\begin{array}{l} g_{x} \\ g_{y} \end{array}\right]\left[\begin{array}{l} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{array}\right]∇f≡grad(f)≡[gxgy][∂x∂f∂y∂f]
∇f\nabla f∇f 表征了 fff 在该点 (x,y)(x,y)(x,y) 处最大变化率方向其大小用 MMM 表示方向角用 α\alphaα 表示方向角是与以左上角为原点的坐标轴 x 轴之间的夹角
M(x,y)gx2gy2M(x, y)\sqrt{g_{x}^{2}g_{y}^{2}}M(x,y)gx2gy2
α(x,y)arctan[gygx]\alpha(x, y)\arctan \left[\frac{g_{y}}{g_{x}}\right]α(x,y)arctan[gxgy]
M(x,y)M(x,y)M(x,y)、α(x,y)\alpha (x,y)α(x,y) 和 fff 的尺寸相同∇f\nabla f∇f 和该点处的边缘正交故梯度向量也称边缘法线。
通过 M(x,y)M(x,y)M(x,y) 和 α(x,y)\alpha (x,y)α(x,y)可以分别构造梯度幅值矩阵和方向矩阵。梯度幅值矩阵是梯度的大小表征了边缘的强度。而梯度方向矩阵则表征了梯度的方向。但是在计算 M(x,y)M(x,y)M(x,y) 时涉及平方和平方根操作计算开销大所以经常使用绝对值来近似梯度幅值
M(x,y)≈∣gx∣∣gy∣M(x, y) \approx\left|g_{x}\right|\left|g_{y}\right|M(x,y)≈∣gx∣∣gy∣
Sobel算子
梯度中偏导数 gxg_xgx、gyg_ygy 的计算需要使用模版算子对 fff 进行滤波。
一维模版对x方向、y方向分别求偏导数即求x方向、y方向的边缘但未考虑对角线方向的偏导数。Roberts[1965]交叉梯度算子以求对角像素45∘,−45∘45^\circ,-45^\circ45∘,−45∘之差为基础但 2×22\times 22×2 大小的模版对于用关于中心点对称的模版来计算边缘方向不是很有用为了考虑中心点对称数据的性质、并携带有关边缘方向的更多信息Prewitt[1970] 算子被提出
Prewittx[−1−1−1000111],Prewitty[−101−101−101]Prewitt_x \begin{bmatrix} -1 -1 -1 \\ 0 0 0 \\ 1 1 1 \end{bmatrix}, Prewitt_y \begin{bmatrix} -1 0 1 \\ -1 0 1 \\ -1 0 1 \end{bmatrix}Prewittx⎣⎡−101−101−101⎦⎤,Prewitty⎣⎡−1−1−1000111⎦⎤
索贝尔算子Sobel operator[1970]主要用作边缘检测本质上是一阶离散性差分算子用来计算图像亮度函数的偏导数的近似值因为它可以反应出灰度值剧烈变化的区域而平滑灰度值变化缓慢的区域。Sobel 算子在中心位置使用 2 以平滑图像更好的抑制噪声Sobel 算子如下
Sobelx[−1−2−1000121],Sobely[−101−202−101]Sobel_x \begin{bmatrix} -1 -2 -1 \\ 0 0 0 \\ 1 2 1 \end{bmatrix}, Sobel_y \begin{bmatrix} -1 0 1 \\ -2 0 2 \\ -1 0 1 \end{bmatrix}Sobelx⎣⎡−101−202−101⎦⎤,Sobely⎣⎡−1−2−1000121⎦⎤
那么
GxSobelx∗fs,GySobely∗fsG_xSobel_x\ast f_s\ ,\ G_ySobel_y\ast f_sGxSobelx∗fs , GySobely∗fs
MGx2Gy2M\sqrt{G_{x}^{2}G_{y}^{2}}MGx2Gy2
αarctan[GyGx]\alpha\arctan \left[\frac{G_{y}}{G_{x}}\right]αarctan[GxGy]
梯度阈值处理
为了减小噪声对梯度幅值矩阵的影响即可以在计算梯度前对图像进行低通滤波平滑也可以在计算梯度后对梯度图像进行阈值处理可以设置大于或等于梯度图像的最大值的 xx% 的像素显示为白色而低于该阈值的像素则显示为黑色以减小噪声的影响。但是在阈值处理后可能会导致部分边缘的断线。
Canny边缘检测核心思想
在Canny[1968]之前Marr和Hildreth[1980]就用高斯拉普拉斯LoG函数生成模版来进行边缘检测他们论证了滤波器 ∇2G\nabla ^2G∇2G 的有效性。
Canny的三个目标
低错误率。所有边缘都应被找到并且没有伪响应边缘点被很好地定位。由检测器标记为边缘的点和真实边缘的中心之间的距离尽可能小单一的边缘点响应。对于真实的边缘点检测器仅返回一个点对于仅存一个单一边缘点的位置检测器不应指出多个边缘像素
Canny基本步骤
高斯滤波平滑输入图像计算梯度幅值图像和角度图像对梯度幅值图像应用非最大抑制NMS双阈值处理划分强、弱边缘滞后边界跟踪消除干扰孤立弱边缘
NMS
非最大值抑制NMS可以细化得到的边缘梯度幅值该方法本质上将边缘按梯度方向的范围进行离散化定义了四个离散方向分别是水平、垂直、45∘45^\circ45∘、−45∘-45^\circ−45∘. 以梯度幅值矩阵上每个点为中心寻找其对于梯度方向最靠近的离散方向若 M(x,y)M(x,y)M(x,y) 至少小于沿该方向的两个邻居之一则抑制gN(x,y)0g_N(x,y)0gN(x,y)0否则 gN(x,y)M(x,y)g_N(x,y)M(x,y)gN(x,y)M(x,y).
NMS 的本质是对各点的由边缘法线确定的边缘进行细化在该点处边缘的两侧和中心选取最大值其余相邻像素位置零。
注意查阅资料得知NMS 也可以使用插值方法来精确计算这里就用本书中提到的方法进行近似计算。 双阈值处理
最后对 gNg_NgN 进行双阈值处理以减少伪边缘点噪声。
在用梯度幅值表征边缘时用单阈值处理梯度幅值若阈值过低则会存在一些伪边缘若阈值过高则会删除有效边缘点造成断线。
改进方式是设置两个阈值LTL_TLT 和 HTH_THT.
创建高、低阈值图像
gNH(x,y)gN(x,y)≥THg_{N H}(x, y)g_{N}(x, y) \ge T_{H}gNH(x,y)gN(x,y)≥TH
gNL(x,y)gN(x,y)≥TLg_{N L}(x, y)g_{N}(x, y) \geq T_{L}gNL(x,y)gN(x,y)≥TL
此时gNHg_{N H}gNH 中的所有非零像素都包含在 gNLg_{N L}gNL 中则做如下运算得到“纯净的” gNLg_{N L}gNL其像素点的灰度值在两阈值之间
gNL(x,y)gNL(x,y)−gNH(x,y)g_{N L}(x, y)g_{N L}(x, y)-g_{N H}(x, y)gNL(x,y)gNL(x,y)−gNH(x,y)
gNHg_{N H}gNH 中的所有强像素均被认定为有效边缘像素立即标记。而由于 THT_HTH 通常较大故 gNHg_{N H}gNH 中的边缘会有缝隙。
滞后边界跟踪
弱边缘点则可能是真的边缘也可能是噪声或颜色变化引起的。为得到精确的结果后者引起的弱边缘点应该去掉。通常认为真实边缘引起的弱边缘点和强边缘点是连通的而噪声引起的弱边缘点则不会。
滞后边界跟踪算法检查一个弱边缘点的 8 连通领域像素只要有强边缘点存在那么这个弱边缘点被认为是真是边缘保留下来。
这个算法搜索所有连通的弱边缘如果一条连通的弱边缘的任何一个点和强边缘点连通则保留这条弱边缘否则抑制这条弱边缘。搜索时可以用广度优先或者深度优先算法。
算法说明如下
输入图像计算其 gNL,gNHg_{NL},g_{NH}gNL,gNH建立visited矩阵和输入图像同形表征每个像素是否在 DFS 中被访问初始化该矩阵为零矩阵建立 gnl_filtered矩阵表征 gNLg_{NL}gNL 中连通的弱边缘且该弱边缘也和 gNHg_{NH}gNH 中任意强边缘连通即过滤调不符合条件的弱边缘后所剩下的符合条件的弱边缘被认为是真边缘。初始化该矩阵为零矩阵不符合条件的弱边缘被抑制置零从上到下自左向右遍历所有强边缘 gNHg_{NH}gNH若该像素已经被访问则继续访问访问下一个像素直到所有强边缘都被访问。若像素未被访问进入 DFS.DFS 中搜索该强边缘的连通弱边缘。先标记该点已被访问如果该点的 gNL(x,y)0g_{NL}(x,y)\gt 0gNL(x,y)0 即该点是弱边缘将gnl_filtered对应位置改为该点的灰度值弱边缘符合条件保留再一次递归搜索该点的8领域寻找连通弱边缘。DFS 的退出条件为已被标记即visited矩阵对应位置为True或访问的点超出边界。输出经过双阈值处理的强边缘矩阵、弱边缘矩阵滞后边界跟踪的连通弱边缘矩阵
之后可对该灰度图像进行二值化并展示。
编码实现和结果分析
import cv2
import numpy as np
import matplotlib.pyplot as pltimg0 house.tif
f cv2.imread(img0, cv2.IMREAD_GRAYSCALE)
plt.imshow(f, cmapgray)
plt.show()def gauss_seperate_kernel(kernel_size, sigma1):create gaussian kernel use separable vector, because gaussian kernel is symmetric and separable.param: kernel_size: input, the kernel size you want to create, normally is uneven number, such as 1, 3, 5, 7param: sigma: input, sigma of the gaussian function, default is 1return a desired size 2D gaussian kernelm kernel_size[0]n kernel_size[1]x_m np.arange(m) - m // 2y_m np.exp(-(x_m ** 2) / (2 * sigma ** 2)).reshape(-1, 1)x_n np.arange(n) - n // 2y_n np.exp(-(x_n ** 2) / (2 * sigma ** 2)).reshape(-1, 1)kernel y_m * y_n.Treturn kerneln 3 # gauss kernel size
sigma 1 # sigma of the gaussian functiongauss_kernel gauss_seperate_kernel((n, n), sigmasigma)
fs cv2.filter2D(f, ddepthcv2.CV_32FC1, kernelgauss_kernel, borderTypecv2.BORDER_DEFAULT)
plt.imshow(fs, cmapgray)
plt.show()sobelx [[-1, -2, -1], [0, 0, 0], [1, 2, 1]]
sobelx np.asarray(sobelx)
sobely sobelx.TGx cv2.filter2D(fs, -1, kernelsobelx, borderTypecv2.BORDER_DEFAULT)
Gy cv2.filter2D(fs, -1, kernelsobely, borderTypecv2.BORDER_DEFAULT)fig, axs plt.subplots(1, 2, figsize(10, 5))
axs[0].set_title(Gx)
axs[0].imshow(Gx, cmapgray)
axs[1].set_title(Gy)
axs[1].imshow(Gy, cmapgray)plt.show()很明显GxG_xGx 表征 xxx 方向上的边界GyG_yGy 表征 yyy 方向上的边界.
M np.sqrt(np.square(Gx) np.square(Gy))
alpha np.arctan(np.divide(Gy, Gx, whereGx!0))fig, axs plt.subplots(1, 2, figsize(10, 5))
axs[0].set_title(M)
axs[0].imshow(M, cmapgray)
axs[1].set_title(r\alpha)
axs[1].imshow(alpha, cmapgray)plt.show()左图为梯度幅值右图为梯度方向。
梯度幅值表征了边缘强度梯度方向较好地识别了背景和房子背景的梯度方向一致。
t .33 # single thresholdMt M.copy()
Mt[Mt t * np.max(Mt)] np.min(Mt)plt.imshow(Mt, cmapgray)
plt.show()可以看到用单阈值直接处理梯度赋值图像得到的边缘中一部分砖瓦的缝隙我们不希望的消失了但也造成了部分边缘断线且此时边缘仍然较粗。
def NMS(gradients, direction): Non-maxima suppressionArgs:gradients: the gradients of each pixeldirection: the direction of the gradients of each pixelReturns:the output imageW, H gradients.shapenms np.copy(gradients)for i in range(1, W - 1):for j in range(1, H - 1):theta direction[i, j]if np.abs(theta) np.pi / 8:d [0, 1] # d1elif np.abs(theta) np.pi * 3 / 8:d [1, 0] # d3elif np.pi / 8 theta:d [1, 1] # d2else:d [1, -1] # d4g1 gradients[i d[0], j d[1]]g2 gradients[i - d[0], j - d[1]]if g1 gradients[i, j] or g2 gradients[i, j]:nms[i, j] 0return nmsMn NMS(M, alpha)# binarization
Mn_b np.copy(Mn)
threshold 0.3 * np.max(Mn_b)
Mn_b[Mn_b threshold] 0
Mn_b[Mn_b threshold] 1fig, axs plt.subplots(1, 2, figsize(10, 5))
axs[0].set_title(Mn)
axs[0].imshow(Mn, cmapgray)
axs[1].set_title(Mn_binarized)
axs[1].imshow(Mn_b, cmapgray)
plt.show()上图为经过 NMS 后的边缘梯度图像边缘被细化了。接下来进行双阈值处理和滞后边缘跟踪进一步细化边缘。
def double_threshold(gn, tl, th): Double ThresholdUse two thresholds to compute the edge.Args:gn: the input imagetl: the low thresholdth: the high thresholdReturns:gnh, gnl, and true_gnlgnh, gnl np.copy(gn), np.copy(gn)tl np.max(gn) * tlth np.max(gn) * thgnh[gnh th] .0gnl[gnl tl] .0gnl gnl - gnhvisited np.zeros_like(gn)gnl_filtered np.zeros_like(gn).astype(float)W, H gn.shapedef dfs(i, j):if i W or i 0 or j H or j 0 or visited[i, j] 1:returnvisited[i, j] 1if gnl[i, j] 0 or gnh[i, j] 0:gnl_filtered[i, j] gnl[i, j]dfs(i - 1, j - 1)dfs(i - 1, j)dfs(i - 1, j 1)dfs(i, j - 1)dfs(i, j 1)dfs(i 1, j - 1)dfs(i 1, j)dfs(i 1, j 1)for w in range(W):for h in range(H):if visited[w, h] 1:continueif gnh[w, h] 0:dfs(w, h)return gnh, gnl, gnl_filtereddef draw_double_threshold(gn, tl, th):draw result pictures:return: gNH_bgNL_true_bgNH, gNL, gNL_true double_threshold(gngn, tltl, thth)gNH_b gNH.copy()gNL_b gNL.copy()gNL_true_b gNL_true.copy()gNH_b[gNH_b 0] np.max(gNH_b)gNL_b[gNL_b 0] np.max(gNL_b)gNL_true_b[gNL_true_b 0] np.max(gNL_true_b)combination gNH_b gNL_true_bcombination[combination 0] np.max(combination)fig, axs plt.subplots(2, 2, figsize(10, 10))axs[0][0].set_title(gNH_b)axs[0][0].imshow(gNH_b, cmapgray)axs[0][1].set_title(gNL_b)axs[0][1].imshow(gNL_b, cmapgray)axs[1][0].set_title(gNL_true_b)axs[1][0].imshow(gNL_true_b, cmapgray)axs[1][1].set_title(combination: gNH_bgNL_true_b)axs[1][1].imshow(combination, cmapgray)plt.suptitle(tl{}, th{}.format(tl, th))plt.show()return combinationtl0 0.05
th0 0.4
binarized_res0 draw_double_threshold(Mn, tl0, th0)tl1 0.1
th1 0.5
binarized_res1 draw_double_threshold(Mn, tl1, th1)适当提高双阈值后(0.05, 0.4)-(0.1, 0.5)对噪声的过滤增强对强边缘的要求提高故强边缘像素减少、弱边缘中噪声也减少。
过滤弱边缘中非连通边缘再将过滤后的弱边缘和强边缘二值化图相加得到最终的二值化边缘结果gNH_bgNL_true_b.
# 和cv2库Canny函数对比
tl2 255 * 0.1
th2 255 * 0.5
g cv2.Canny(f, tl2, th2)
plt.imshow(g, cmapgray)
plt.show()参考文献
数字图像处理第3版北京电子工业出版社图像滤波之高斯滤波边缘检测之Sobel检测算子Canny边缘检测算法的python实现常用测试数据集下载“Standard” test images
