图像边缘检测(canny检测附完整代码与输出样例)_图像边缘检测代码

首先通过用高斯核对输入图像进行卷积来平滑输入图像。

大小为(2k+1) × (2k+1)的高斯核的方程由下式给出：

导入所需模块：

from __future__ import print_function
 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from time import time
from skimage import io
 
plt.rcParams['figure.figsize'] = (15.0, 12.0) # 设置绘图默认大小
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'
 
%load_ext autoreload #自动加载拓展板块
%autoreload 2

定义高斯核：

import numpy as np
def gaussian_kernel(size, sigma):
    kernel = np.zeros((size, size))
    if size%2 == 0:
        print("not odd")
    else:
        k=(size-1)/2
        for i in range(size):
            for j in range(size):
                kernel[i,j] = np.exp(-((i-k)**2 + (j-k)**2)/(2*sigma**2))/(2*np.pi*sigma**2)
    pass
    return kernal

定义卷积运算：

def conv(image, kernel):
    Hi, Wi = image.shape    #图像尺寸
    Hk, Wk = kernel.shape    #卷积核大小
    out = np.zeros((Hi, Wi))    #输出尺寸
    pad_width0 = Hk // 2    #纵向填充
    pad_width1 = Wk // 2    #横向填充
    pad_width = ((pad_width0,pad_width0),(pad_width1,pad_width1))
    padded = np.pad(image, pad_width, mode='edge')    #填充操作
    
    kernel=np.flip(kernel)   #卷积会翻转内核，此时将卷积核翻转
    for i in range(Hi):
        for j in range(Wi):
            temp=padded[i:Hk+i,j:Wk+j]    #填充操作
            out[i,j]=np.sum(temp*kernel)    #计算每个像素的邻域加权和
    pass
    return out

任选一张图片（doge）作为测试样例

对图像平滑处理

kernel_size = 5 
sigma = 1.4    #可改变参数对比图像模糊程度
 
img = io.imread('dog.png', as_gray=True)    #加载需处理的图像
kernel = gaussian_kernel(kernel_size, sigma)    #使用高斯核作为卷积核
smoothed = conv(img, kernel)    #对图像平滑处理
#绘图操作
plt.subplot(1,2,1)
plt.imshow(img)
plt.title('Original image')
plt.axis('off')
 
plt.subplot(1,2,2)
plt.imshow(smoothed)
plt.title('Smoothed image')
plt.axis('off')
 
plt.show()

寻找梯度：

实现x 与y 方向上的导数：

def partial_x(img):
    out = None
    kernel=np.array([[1,0,-1]])/2    #定义卷积核
    out=conv(img,kernel)
    pass
    return out
 
def partial_y(img):
    out = None
    kernel=np.array([[1],[0],[-1]])/2
    out=conv(img,kernel)
    pass
    return out

计算平滑处理后图像在x，y方向上的偏导数

Gx = partial_x(smoothed)
Gy = partial_y(smoothed)
 
plt.subplot(1,2,1)
plt.imshow(Gx)
plt.title('Derivative in x direction')
plt.axis('off')
 
plt.subplot(1,2,2)
plt.imshow(Gy)
plt.title('Derivative in y direction')
plt.axis('off')
 
plt.show()

输出在x与y方向上导数图像

求图像梯度大小与方向

利用上述计算的偏导数可以计算梯度的大小和方向，其公式如下所示

def gradient(img):
    G = np.zeros(img.shape)
    theta = np.zeros(img.shape)
    Gx=partial_x(img)
    Gy=partial_y(img)
    G=np.sqrt(Gx**2+Gy**2)    #计算梯度的大小
    theta=np.arctan2(Gy,Gx)    #计算梯度的方向
    theta=theta*180/np.pi%360    #将方向转化为角度
    pass
    return G, theta

G, theta = gradient(smoothed)
 
if not np.all(G >= 0):    #梯度大小应为非负数
    print('Magnitude of gradients should be non-negative.')
    
if not np.all((theta >= 0) * (theta < 360)):    #梯度方向应在0<=θ<360的范围内
    print('Direction of gradients should be in range 0 <= theta < 360')
 
plt.imshow(G)
plt.title('Gradient magnitude')
plt.axis('off')
plt.show()

输出梯度图像

非极大值抑止

此步骤的目的是将“模糊”边缘转换为“锐利”边缘。基本上，这是通过保留梯度图像中的所有局部最大值并丢弃其他所有内容来完成的。该算法针对梯度图像中的每个像素 (x,y)：

1. 将梯度方向四舍五入到最近的45度，对应使用8-连通邻域

2.将当前像素的边缘强度与该像素在正负梯度方向的边缘强度进行比较。例如，如果梯度方向为南（theta=90），则与南北方向的像素进行比较。

3如果当前像素的边缘强度最大；保留边缘强度的值。如果不是则删除该值。

def non_maximum_suppression(G, theta):
    H, W = G.shape
    out = np.zeros((H, W))
 
    theta = np.floor((theta + 22.5) / 45) * 45    #将渐变方向四舍五入到最接近的45度
    G = np.pad(G,((1,1),(1,1)),mode = 'constant')
 
    for i in range(1,H+1):
        for j in range(1,W+1):
            # 计算方向
            rad = np.deg2rad(theta[i-1, j-1])#角度转弧度计算sincos
            h =int(np.around(np.sin(rad)))   # 行
            w =int(np.around(np.cos(rad)))   # 列
            #得到当前位置以及梯度方向的两个相邻值
            p1 = G[i+h, j+w]
            p2 = G[i-h, j-w]
            if(G[i, j] > p1 and G[i, j] > p2): # 如果当前位置是该方向的最大值，则保留
                out[i-1, j-1] = G[i, j]
            else:
                out[i-1, j-1] = 0
    pass
    return out

输出非极大抑制图像

nms = non_maximum_suppression(G, theta)
plt.imshow(nms)
plt.title('Non-maximum suppressed')
plt.axis('off')
plt.show()

可以看出该图像的边缘与梯度图像相比锐利程度更高。

双阈值处理

在非最大抑制步骤之后剩余的边缘像素仍然用它们的强度逐像素标记。 其中许多可能是图像中的真实边缘，但有些可能是由噪声或颜色变化引起的。区分这些之间的最简单方法是使用阈值，以便仅保留某个值的强边缘。 Canny 边缘检测算法使用双阈值。 比高阈值强的边缘像素被标记为强； 弱于低阈值的边缘像素被抑制，两个阈值之间的边缘像素被标记为弱。

定义双阈值算法：

def double_thresholding(img, high, low):
    strong_edges = np.zeros(img.shape, dtype=np.bool)
    weak_edges = np.zeros(img.shape, dtype=np.bool)
    H,W=img.shape
    for i in range(H):
        for j in range(W):
            if img[i,j]>high:    #像素大于阈值时
                strong_edges[i,j]=img[i,j]    #强边缘被赋值
            elif img[i,j]<high and img[i,j]>low:    #像素小于阈值时
                weak_edges[i,j]=img[i,j]    #弱边缘被赋值
    pass
    return strong_edges,weak_edges

我们自己设置强弱阈值，并对输出图像进行对比

low_threshold = 0.025    #设置低阈值
high_threshold = 0.05    #设置高阈值
 
strong_edges, weak_edges = double_thresholding(nms, high_threshold, low_threshold)
assert(np.sum(strong_edges & weak_edges) == 0)
 
edges=strong_edges * 1.0 + weak_edges * 0.5
 
plt.subplot(1,2,1)
plt.imshow(strong_edges)
plt.title('Strong Edges')
plt.axis('off')
 
plt.subplot(1,2,2)
plt.imshow(edges)
plt.title('Strong+Weak Edges')
plt.axis('off')
 
plt.show()

边缘追踪

强边缘被解释为“某些边缘”，可以立即包含在最终的边缘图像中。如果它直接通过低和高之间的像素连接到“强边缘像素”，则迭代地考虑其相邻像素然后将其声明为“边缘像素”。因为噪声和其他小的变化不太可能导致强边缘（适当调整阈值水平）。因此，强边缘（几乎）仅归因于原始图像中的真实边缘。弱边缘可能是由于真实边缘或噪声/颜色变化。后一种类型可能会在整个图像上独立于边缘分布，因此只有少量位于强边缘附近。由于真实边缘而导致的弱边缘更有可能直接连接到强边缘

我们首先要实现对相邻像素之间的追踪

def get_neighbors(y, x, H, W):
    neighbors = []
    for i in (y-1, y, y+1):
        for j in (x-1, x, x+1):
            if i >= 0 and i < H and j >= 0 and j < W:    #确保迭代次数不超出图像边界
                if (i == y and j == x):
                    continue
                neighbors.append((i, j))
 
    return neighbors

在strong_edges中的每个像素上迭代，并在weak_edges的连接像素上执行广度优先搜索，以将它们链接起来。

def link_edges(strong_edges, weak_edges):
    H, W = strong_edges.shape
    indices = np.stack(np.nonzero(strong_edges)).T
    edges = np.zeros((H, W), dtype=np.bool)x
 
    weak_edges = np.copy(weak_edges)
    edges = np.copy(strong_edges)    #创建新的实例，保存之前的参数
 
    visited = np.zeros((H, W))
    while len(indices) != 0:#强边缘没被删光，这个循环类似于栈，不断把强边缘的非零值以及未被访问过的和强边缘连通的弱边缘赋给edges
        i, j = indices[0]
        edges[i][j] = True
        indices = np.delete(indices, 0, axis=0)
        neighbors = get_neighbors(i, j, H, W)    #追踪相邻像素
        for neighbor in neighbors:
            if weak_edges[neighbor] and not visited[neighbor]:#弱边缘和强边缘连通且未被访问
                visited[neighbor] = True
                indices = np.vstack((indices, neighbor))
    return edges

实现边缘追踪

edges = link_edges(strong_edges, weak_edges)
 
plt.imshow(edges)
plt.axis('off')
plt.show()

canny边缘检测器

将上述的图像处理方法综合为canny边缘检测器

def canny(img, kernel_size, sigma, high, low):
    kernel = gaussian_kernel(kernel_size, sigma)    #选择卷积核为高斯核
    smoothed = conv(img, kernel)    #对图像平滑处理
    G, theta = gradient(smoothed)    #求梯度进行非极大值抑制
    nms = non_maximum_suppression(G, theta)    #非极大值抑制
    strong_edges, weak_edges = double_thresholding(nms, high, low)#双阈值处理
    edge = link_edges(strong_edges, weak_edges)    #边缘追踪
    pass
    return edge