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

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

导入所需模块：

from __future__ import print_functionimport numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from time import time
from skimage import ioplt.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)/2for 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)passreturn 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)    #计算每个像素的邻域加权和passreturn 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 = Nonekernel=np.array([[1,0,-1]])/2    #定义卷积核out=conv(img,kernel)passreturn outdef partial_y(img):out = Nonekernel=np.array([[1],[0],[-1]])/2out=conv(img,kernel)passreturn 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    #将方向转化为角度passreturn 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.shapeout = 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])#角度转弧度计算sincosh =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] = 0passreturn 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.shapefor 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]    #弱边缘被赋值passreturn 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.5plt.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):continueneighbors.append((i, j))return neighbors

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

def link_edges(strong_edges, weak_edges):H, W = strong_edges.shapeindices = np.stack(np.nonzero(strong_edges)).Tedges = np.zeros((H, W), dtype=np.bool)xweak_edges = np.copy(weak_edges)edges = np.copy(strong_edges)    #创建新的实例，保存之前的参数visited = np.zeros((H, W))while len(indices) != 0:#强边缘没被删光，这个循环类似于栈，不断把强边缘的非零值以及未被访问过的和强边缘连通的弱边缘赋给edgesi, j = indices[0]edges[i][j] = Trueindices = 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] = Trueindices = 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)    #边缘追踪passreturn edge