Python图像处理【18】边缘检测详解

0. 前言

边缘是图像中两个区域之间具有相对不同灰级特性的边界，或者说是亮度函数突然变化的像素。边缘检测器是一组重要的局部图像预处理方法，用于定位强度函数中的(急剧)变化。可以通过使用以下算法来检测边缘上的点：

• 一阶导数的局部最大值或最小值
• 二阶导数的零交叉点 (zero-crossing)

判断一个边缘检测算法是否足够优秀可以使用以下标准：

• 能够最大程度地减少假阳性(即虚假边缘)和假阴性(即缺失真实边缘)的概率
• 检测到的边缘必须尽可能接近真正边缘

在本节中，首先，我们将学习如何使用二阶导数的过零点检测图像中的边缘。然后，介绍如何使用深度学习模型进行边缘检测，并将其与 Canny 边缘检测器进行比较。

1. 图像导数

我们首先介绍如何计算图像导数（梯度）以及边缘对图像导数的影响（以便我们可以使用图像导数进行边缘检测）。我们可以使用特定核/掩码的卷积操作近似图像导数：

• 图像一阶导，掩码 $\left[\begin{array}{cc}?1& 1\end{array}\right]$
  $$f^{\prime }(x)=\underset{x\to \infty }{\lim }\frac{f(x+h)?f(x)}{h}=f(x+1)?f(x)$$
• 图像二阶导，掩码$\left[\begin{array}{ccc}1& ?2& 1\end{array}\right]$
  $$\begin{array}{rl}{f}^{\prime \prime }(x)& =\underset{x\to \infty }{\lim }\frac{f^{\prime }(x+h)?f^{\prime }(x)}{h}=f^{\prime }(x)?f^{\prime }(x?1)\\ & =f(x+1)?2f(x)+f(x?1)\end{array}$$

使用以上近似函数计算二值图像的导数，并注意其对导数的影响，以便理解为什么导数可以用于边缘检测。

(1) 首先导入所需的模块、函数，然后读取二值黑白图像：

from scipy.signal import convolve
from skimage.io import imread
from skimage.color import rgb2gray
import matplotlib.pylab as plt

img = rgb2gray(imread('example.png'))
h, w = img.shape

(2) 构建卷积核，使用空间卷积计算一阶导数 $\frac{?f}{?x}$ 和二阶导数 $\frac{{?}^{2}f}{{?}^{2}x}$：

kd1 = [[1, -1]]
kd2 = [[1, -2, 1]]
imgd1 = convolve(img, kd1, mode='same')
imgd2 = convolve(img, kd2, mode='same')

(3) 最后，绘制输入图像及其导数：

plt.figure(figsize=(20,10))
plt.gray()
plt.subplot(231), plt.imshow(img), plt.title('image', size=10)
plt.subplot(232), plt.imshow(imgd1), plt.title('1st derivative', size=10)
plt.subplot(233), plt.imshow(imgd2), plt.title('2nd derivative', size=10)
plt.subplot(234), plt.plot(range(w), img[0,:]), plt.title('image function', size=10)
plt.subplot(235), plt.plot(range(w), imgd1[0,:]), plt.title('1st derivative function', size=10)
plt.subplot(236), plt.plot(range(w), imgd2[0,:]), plt.title('2nd derivative function', size=10)
plt.show()

从以上结果图像可以看出，在边缘像素周围可以观察到以下内容：

• 原始图像中的边缘像素强度急剧变化
• 一阶导数在边缘像素处达到最大值
• 二阶导数在边缘像素处具有零交叉点

2. LoG/zero-crossing

在本节中，我们将使用导数的零交叉 (zero-crossing) 特性查找图像中的边缘。在边缘像素上，一阶导数达到最大化(或最小化)，而此处的二阶导数为零。然而，我们并不能总是找到导数为零的离散像素，因此需要寻找零交叉，以近似与梯度最大值/最小值相对应的位置。
但是，这种方法对噪声较为敏感(因为它需要两次求导)，为了解决这个问题，我们需要首先平滑图像并去除噪声。因此，有以下两种方法使用二阶导数来识别边缘：

• 首先执行平滑，然后应用梯度
• 结合平滑和梯度操作

导数是使用拉普拉斯 (Laplacian $?2$) 算子计算的，并使用高斯算子平滑图像。这两个算子通常可以被组合为高斯拉普拉斯算子 (Laplacian of Gaussian, LoG)，因此可以减少卷积运算：
卷积的导数定理可以描述为：

$$\frac{?}{?x}(h?f)=(\frac{?}{?x}h)?f$$

高斯拉普拉斯算子 (Laplacian of Gaussian, LoG)，或称 Marr-Hildteth 算子描述如下：

$${?}^2[f(x,y)?G(x,y)]={?}^{2}G(x,y)?f(x,y)$$

其中 $G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{?\frac{x^2+y^2}{2{\sigma }^2}}$，将 $G(x,y)$ 带入上式可得：

$$\begin{array}{rl}LoG(x,y)& ={?}^{2}G(x,y)=\frac{{?}^2}{?x^2}G(x,y)+\frac{{?}^2}{?y^2}G(x,y)\\ & =?\frac{1}{\pi {\sigma }^4}[1?\frac{x^2+y^2}{2{\sigma }^2}]e^{?\frac{x^2+y^2}{2{\sigma }^2}}\end{array}$$

参数 $\sigma$ 是高斯核的宽度，并用于控制平滑量。紧接着图像的 LoG 变换之后，对零交叉的连续计算会输出图像中的边缘，这通常称为 Marr-Hildteth (LoG) 算法。

2.1 Marr-Hildteth 算法

Marr-Hildteth (LoG) 算法描述如下：

• 计算图像 LoG
• 在每一行和每一列中找到零交叉点
• 计算零交叉点的斜率
• 将阈值应用于梯度并标记边缘

使用该算法检测到的边缘结构根据高斯核宽度参数 $\sigma$ 不同而不同：

• 检测到更大尺度的边缘，减少噪声，但增加边缘位置的不确定性
• 检测到更精细的细节特征

但该算法不能很好地处理角点；LoG 的零交叉点能够对边缘进行较好的定位，尤其是当边缘不是很锐利时，该方法存在噪声抑制。

(1) 导入所有必需的库，读取输入图像，并根据 LoG/zero-crossings 找到其中的边缘：

import numpy as np
from scipy import ndimage
from skimage.io import imread
import matplotlib.pyplot as plt
from skimage.color import rgb2gray

(2) 定义函数 any_neighbor_neg()，该函数根据给定(正值)像素，返回其周围 8 个邻居中是否有负值像素，如果存在负值像素，则意味着存在零交叉点，将计算并返回相应边的斜率：

def any_neighbor_neg(img, i, j):
    for k in range(-1,2):
        for l in range(-1,2):
            if img[i+k, j+k] < 0:
                return True, img[i, j] - img[i+k, j+k]
    return False, None

(3) 定义函数 zero_crossing()，该函数以黑色输出图像开始，如果以下条件均为真，则将像素颜色更改为白色，即边缘像素：

• 像素是一个正像素
• 它的 8 个邻居中至少有一个负像素
• 边缘的斜率(如果存在负像素)高于给定的阈值(仅检测强边缘)

def zero_crossing(img, th):
    out_img = np.zeros(img.shape)
    for i in range(1,img.shape[0]-1):
        for j in range(1,img.shape[1]-1):
            found, slope = any_neighbor_neg(img, i, j)
            if img[i,j] > 0 and found and slope > th:
                out_img[i,j] = 255
    return out_img

(4) 读取输入图像并将其转换为灰度图像，使用不同 $\sigma$ 值调用 scipy.ndimage 模块的 gaussian_lapace() 函数将 LoG 运算符应用于图像。使用不同 $\sigma$ 值以及不同阈值计算所得图像中的零交叉点，随着 $\sigma$ 的增加，使用的阈值会降低，并绘制使用不同 $\sigma$ 检测到到的的边缘：

img = rgb2gray(imread('1.png'))
#img = misc.imread('../new images/tagore.png')[...,3]
print(np.max(img))
fig = plt.figure(figsize=(10,16))
plt.subplots_adjust(0,0,1,0.95,0.05,0.05)
plt.gray() # show the filtered result in grayscale
for sigma, thres in zip(range(3,10,2), [1e-3, 1e-4, 1e-5, 1e-6]):
    plt.subplot(2,2,sigma//2)
    result = ndimage.gaussian_laplace(img, sigma=sigma)
    result = zero_crossing(result, thres)
    plt.imshow(result)
    plt.axis('off')
    plt.title('LoG with zero-crossing, sigma=' + str(sigma), size=10)

plt.tight_layout()
plt.show()

3. Canny 与 holistically-nested 算法

在本节中，我们将学习一种基于深度学习的边缘检测技术，称为 holistically-nested 边缘检测。我们还将结果与另一种主流的基于梯度的边缘检测算法 Canny 进行比较，我们首先将介绍这两种不同的算法。

3.1 Canny 边缘检测

Canny 边缘检测是一种多阶段算法，依赖于找到图像导数的极值。接下来，我们将使用 OpenCV 实现 Canny 边缘检测：

• 降噪：作为一种基于梯度的边缘检测技术，它易受图像中的噪声的影响；因此，第一步是用高斯滤波器去除图像中的噪声(例如，使用核尺寸 5 x 5)
• 然后，用水平和垂直方向的 Sobel 核(使用卷积)对图像平滑图像的强度梯度计算执行滤波，Sobel 垂直和水平核如下，本质上只是导数算子的变体：
  $$S_x=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ ?1& ?2& ?1\end{array}\right],S_y=\left[\begin{array}{ccc}+1& 0& ?1\\ +2& 0& ?2\\ +1& 0& ?1\end{array}\right]$$
  根据以上公式，每个像素的边缘梯度和方向计算如下，梯度方向始终垂直于边缘：
  • 梯度大小： $|G|=\sqrt{S_x^2+S_y^2}$
  • 角度：$\theta =ta{n}^{?1}\frac{S_x}{S_y}$
• 非极大值抑制：在计算梯度大小和方向之后，对图像进行全扫描以去除可能不构成边缘的任何不需要的像素。为此，在每个像素上，检查该像素在梯度方向上是否是其邻域中的局部最大值，得到具有较小边缘的二值图像
• 滞后性：这个阶段在决定在以上步骤中检测到的边缘中，哪些是真实的边缘，哪些边缘是虚假的边缘。为此，使用阈值 minVal 和 maxVal，强度大于 maxVal 的边缘都是确定的真实边缘， 而小于 minVal 的边缘都是确定的虚假边缘，因此需要将其丢弃。而位于这两个阈值之间的边缘，如果它们连接到确定的边缘像素，则认为它们是边缘的一部分；否则，它们也会被丢弃。

3.2 holistically-nested 边缘检测

Holistically-nested 边缘检测 (Holistically-nested Edge Detection, HED) 是一种新的边缘检测算法，该算法解决了边缘检测算法中的两个重要问题：

• 端到端的模型训练和预测
• 多尺度多层次特征学习

该模型使用深度学习模型执行端到端的预测，可以充分利用全卷积神经网络和深度监督网络。HED 会自动学习丰富的层次表示，它比基于 CNN 的边缘检测算法快数个数量级。
首先下载预处理模型，并保存在 models 文件夹中，网络架构如下所示：

(1) 导入所需库和函数：

import cv2
import numpy as np
import matplotlib.pylab as plt

(2) 加载输入图像并获取其尺寸：

image = cv2.imread('2.png')
(h, w) = image.shape[:2]

(3) 将图像转换为灰度图像，用高斯模糊进行模糊，并使用 Canny 边缘检测获得边缘，滞后阈值分别为 80 和 150：

gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
blurred = cv2.GaussianBlur(gray, (5, 5), 0)
canny = cv2.Canny(blurred, 80, 150) 

(4) 接下来，我们将用一个居中的裁剪层替换 OpenCV 的裁剪层，使用 getMemoryShapes() 和 forward() 方法创建类 CropLayer()。在 getMemoryShapes() 方法中，裁剪层接收两个输入，并保持批大小和通道数量，计算开始和结束裁剪坐标，并返回数据集形状；在 forward() 方法中，使用 divide(x) 执行裁剪：

class CropLayer(object):
    def __init__(self, params, blobs):
        self.xstart = 0
        self.xend = 0
        self.ystart = 0
        self.yend = 0
        
    def getMemoryShapes(self, inputs):
        inputShape, targetShape = inputs[0], inputs[1]
        batchSize, numChannels = inputShape[0], inputShape[1]
        height, width = targetShape[2], targetShape[3]
        self.ystart = (inputShape[2] - targetShape[2]) // 2
        self.xstart = (inputShape[3] - targetShape[3]) // 2
        self.yend = self.ystart + height
        self.xend = self.xstart + width
        return [[batchSize, numChannels, height, width]]
    
    def forward(self, inputs):
        return [inputs[0][:,:,self.ystart:self.yend,self.xstart:self.xend]]

(5) 接下来，读取预训练的模型：

prototxt_path = "models/deploy.prototxt"
model_path = "models/hed_pretrained_bsds.caffemodel"

net = cv2.dnn.readNetFromCaffe(prototxt_path, model_path)

(6) 使用 dnn_registerLayer() 将新层注册到模型中：

cv2.dnn_registerLayer('Crop', CropLayer)

(7) 使用 blobFromImage() 为 HED 模型从输入图像构造输入 blob:

blob = cv2.dnn.blobFromImage(image, scalefactor=1.0, size=(w, h), mean=(104.00698793, 116.66876762, 122.67891434),  swapRB=False, crop=False)

(8) 将 blob 设置为网络的输入，并运行正向传播 format 方法计算输出：

net.setInput(blob)
outs = net.forward()
hed = cv2.resize(outs[0][0,:,:], (w, h))
hed = (255 * hed).astype("uint8")

(9) 通过 Canny 和 Holistically-nested 边缘检测获得的输出边缘检测结果如下所示：

plt.figure(figsize=(20, 8))
plt.gray()
plt.subplots_adjust(0,0,1,0.975,0.05,0.05)
plt.subplot(131), plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB)), plt.axis('off'), plt.title('input', size=10)
plt.subplot(132), plt.imshow(canny), plt.axis('off'), plt.title('canny', size=10)
plt.subplot(133), plt.imshow(hed), plt.axis('off'), plt.title('holistically-nested', size=10)
plt.show()

小结

边缘检测是图像处理和计算机视觉中的基本问题，边缘检测的目的是标识图像中亮度变化明显的点，图像属性中的显著变化通常反映了图像中的重要变化和特征。边缘检测是图像处理和计算机视觉中，尤其是特征提取中的一个重要研究领域。在本节中，我们学习了多种边缘检测算法，包括基于梯度的算法以及基于深度神经网络的方法。

系列链接

Python图像处理【1】图像与视频处理基础
Python图像处理【2】探索Python图像处理库
Python图像处理【3】Python图像处理库应用
Python图像处理【4】图像线性变换
Python图像处理【5】图像扭曲/逆扭曲
Python图像处理【6】通过哈希查找重复和类似的图像
Python图像处理【7】采样、卷积与离散傅里叶变换
Python图像处理【8】使用低通滤波器模糊图像
Python图像处理【9】使用高通滤波器执行边缘检测
Python图像处理【10】基于离散余弦变换的图像压缩
Python图像处理【11】利用反卷积执行图像去模糊
Python图像处理【12】基于小波变换执行图像去噪
Python图像处理【13】使用PIL执行图像降噪
Python图像处理【14】基于非线性滤波器的图像去噪
Python图像处理【15】基于非锐化掩码锐化图像
Python图像处理【16】OpenCV直方图均衡化
Python图像处理【17】指纹增强和细节提取