昨天 10:22 常州大学 C++ 发布于四川

关注

OpenCV边缘检测算法数学原理

边缘是指图像中灰度发生急剧变化的区域。两个具有不同灰度值的相邻区域之间总存在着边缘，这是灰度值不连续的结果。这种不连续通常可以利用求导的方法检测到，一般常用一阶、二阶和一二阶混合导数来检测边缘。边缘检测的数学原理需要一些基础数学知识。

OpenCV内置边缘检测算法

Sobel算子

Sobel算子的原理可参照此篇文章：关于Sobel核的一般数学推导

Scharr算子

Scharr算子卷积核为：

$\begin{aligned} Scharr_X&=\begin{bmatrix} -3&0&3\\ -10&0&10\\ -3&0&3 \end{bmatrix}\\ Scharr_Y&=\begin{bmatrix} -3&-10&-3\\ 0&0&0\\ 3&10&3 \end{bmatrix} \end{aligned}$

Scharr算法在OpenCV中有Sobel接口（将ksize=1）和Scharr接口两种调用方式，为语义清晰，推荐后一种方式。

Laplacian算子

Laplacian边缘检测算法是二阶算法，表达式为

$\begin{aligned} \nabla^2f(x,y)&=\frac{\partial^2f(x,y)}{\partial x^2}+\frac{\partial^2f(x,y)}{\partial y^2}\\ &=f(x+1,y)-2f(x,y)+f(x-1,y)+f(x,y+1)-2(x,y)+f(x,y-1)\\ &=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y) \end{aligned}$

这里需要注意对二阶梯度进行区分。梯度以矩阵表示，Laplacian以标量表示。这样我们得到了Laplacian的第一个卷积核：

$Laplacian=\begin{bmatrix} 0&1&0\\ 1&-4&1\\ 0&1&0 \end{bmatrix}$

因为上式仅考虑了水平和垂直方向的情况，实践中也给出了另一个Laplacian卷积核：

$Laplacian=\begin{bmatrix} 1&1&1\\ 1&-8&1\\ 1&1&1 \end{bmatrix}$

Canny算子

Canny边缘检测算法包含以下五个主要步骤，具有检测率高、定位准确、一次响应三个特点。

噪声抑制

Canny算法第一步是使用高斯滤波来去除图像中的噪声。设原图像为 $f(x,y)$ ，矩阵为 $F$ 。高斯滤波的核矩阵为 $K$ ，则第一步得出图像与高斯滤波的卷积 $F*K$ 。需要注意的是，此步并未包含在OpenCV中的Canny接口中，需要自行调用接口。当然，这也意味着此步可以根据具体条件使用除了高斯滤波之外的其他滤波方法。本文中仅考虑线性滤波的情况，非线性滤波不适用于该表达（读者可自行推导）。

计算梯度

使用Sobel算子 $S_x,S_y$ ，得到 $x$ 方向和 $y$ 方向的 $G_x,G_y$ ，则梯度为

$G=\begin{bmatrix} G_x & G_y \end{bmatrix}=\begin{bmatrix} F*K*S_x & F*K*S_y \end{bmatrix}$

梯度幅值有以下两个表达式，可根据特点来自行选取：

$\begin{aligned} mag(G)&=|G_x|+|G_y|\\ mag(G)&=\sqrt{G_x^2+G_y^2} \end{aligned}$

梯度方向：

$\theta=\arctan\left(\frac{G_y}{G_x}\right)$

非极大值抑制

非极大值抑制(Non-Maximum Suppression, NMS)是使边缘细化为单像素宽度的关键步骤，目的是消除在梯度方向上不是真正边缘的像素，从而得到更加精确的边界定位。NMS步骤：

对每个像素点，沿梯度方向 $\theta$ ，找出前后相邻的两个像素的幅值 $mag(G_1),mag(G_2)$ 和当前像素点幅值 $mag(G)$
如果当前点的幅值不是该方向上的最大值（ $mag(G)\not=max(mag(G),mag(G_1),mag(G_2))$ ），则将其抑制（设为0）

方向角 $\theta$ 取值为 $45^\circ$ 的整数倍，通常近似取 $0^\circ(水平),45^\circ(对角),90^\circ(垂直),135^\circ(对角)$ 。如计算出的梯度方向角 $\theta=30^\circ$ ，由于离 $45^\circ$ 更近，则近似取 $\theta\approx45^\circ$ 。需要注意的是，关于幅度的计算，若

$mag_1(G_1)=|G_{1x}|+|G_{1y}|<mag_1(G_2)=|G_{2x}|+|G_{2y}|$

并不足以说明

$mag_2(G_2)=\sqrt{G_{1x}^2+G_{1y}^2}<mag_2(G_1)=\sqrt{G_{2x}^2+G_{2y}^2}$

可举反例说明： $1+100=101<102=60+60$ ，然而 $\sqrt{1^2+100^2}=\sqrt{10001}>\sqrt{7200}=\sqrt{60^2+60^2}$ 。也就是说，在需要高精度的场景，应当使用平方和的开方，而不是绝对值的和；需要高效的场景反之。

双阈值检测

定义两个阈值 $T_H,T_L$ ，将幅值与两个阈值对比，根据对比情况依次分为强边缘、弱边缘和非边缘。我们保留强边缘，滞后处理弱边缘，抑制非边缘。具体为

强边缘： $mag(G)>T_H$
弱边缘： $T_L<mag(G)\le T_H$
非边缘： $mag(G)\le T_L$

边缘连接

对于上一步的弱边缘，判断其是否与强边缘连接。如果是，则保留；否则抑制。此步从所有强边缘像素点出发，使用广度优先搜索（BFS）算法，将所有相连的弱边缘像素点转换为强边缘像素点。其余未被搜索到的弱边缘点被抑制。

Canny算法对真实边缘进行细化，较好的排除了非边缘的影响（如色彩渐变）。

其他边缘检测算法

差分边缘检测

差分边缘检测是最基础的边缘检测算法。差分边缘检测算法原理是：将图像的相邻像素点进行差分，得到像素点的梯度，然后根据梯度的大小来判断像素点是否为边缘点。差分边缘检测在 $x$ 方向上的一阶差分定义为 $G_x(x,y)=f(x+1,y)-f(x,y)$ ，算子定义为

$Diff_x=\begin{bmatrix} 0&0&0\\ 1&-1&0\\ 0&0&0 \end{bmatrix}$

在 $y$ 方向上的一阶差分定义为 $G_y(x,y)=f(x,y+1)-f(x,y)$ ，核模板为

$Diff_y=\begin{bmatrix} 0&1&0\\ 0&-1&0\\ 0&0&0 \end{bmatrix}$

前者用于增强垂直边边缘，后者增强计算水平边缘。

Roberts算子

Roberts边缘检测算法又称为交叉微分算法，是基于交叉微分的梯度算法。Roberts算子的核模板为

$Roberts_x=\begin{bmatrix} -1&0\\ 0&1 \end{bmatrix} ,Roberts_y=\begin{bmatrix} 0&-1\\ 1&0 \end{bmatrix}$

差分定义为

$\begin{aligned} G_x(x,y)&=f(x+1,y+1)-f(x,y)\\ G_y(x,y)&=f(x,y+1)-f(x+1,y) \end{aligned}$

Roberts算子能较好地增强 $\pm45^\circ$ 的图像边缘。

Prewitt算子

上述两种算子由于每次只关注了两个像素，因此容易受到噪声干扰。Prewitt每次基于边缘的6个像素的平均值来计算梯度，因此可以减少噪声，但同时对边缘的定位不如上两个算法。核模板为

$\begin{aligned} Prewitt_x&=\begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix}\cdot\begin{bmatrix} -1&0&1 \end{bmatrix}=\begin{bmatrix} -1&0&1\\ -1&0&1\\ -1&0&1 \end{bmatrix}\\ Prewitt_y&=\begin{bmatrix} 1\\ 0\\ -1 \end{bmatrix}\cdot\begin{bmatrix} 1&1&1 \end{bmatrix}=\begin{bmatrix} 1&1&1\\ 0&0&0\\ -1&-1&-1 \end{bmatrix} \end{aligned}$

从分离的结果来看， $Prewitt_x$ 算子实际上是先对图像进行垂直方向的非归一化的均值平滑，在进行水平方向的差分； $Prewitt_y$ 算子则先对图像进行水平方向的非归一化的均值平滑，在进行垂直方向的差分。这就是Prewitt算子能够抑制噪声的原因。 Prewitt算子还包含两个对角线版本的算子：

$\begin{bmatrix} 0&1&1\\ -1&0&1\\ -1&-1&0 \end{bmatrix},\begin{bmatrix} -1&-1&0\\ -1&0&1\\ 0&1&1 \end{bmatrix}$

LoG算子

LoG（Laplacian-of-Gaussian）边缘检测算法首先对图像进行高斯滤波，然后求其拉普拉斯二阶导数，即图像与Laplacian of Gaussian方法进行滤波运算。最后通过检测滤波结果的零交叉以获得图像或物体的边缘。

滤波

LoG算法第一步是使用高斯滤波来去除图像中的噪声。设原图像为 $f(x,y)$ ，矩阵为 $F$ 。高斯滤波的核矩阵为 $K$ ，则第一步得出图像与高斯滤波的卷积 $F*K$ 。这里跟Canny算法第一步是相似的，不同之处在于Canny滤波可以使用其他滤波方法，这里由于算法本身的定义，只采用高斯滤波。

增强

对上一步处理后的平滑图像进行Laplacian运算，即 $h(x,y)=\nabla^2((f*k)(x,y))$ 。对平滑图像进行Laplacian运算可等效为先对高斯滤波卷积核进行Laplacian运算，再与 $f(x,y)$ 卷积，即 $h(x,y)=(f*\nabla^2k)(x,y)$ 。其中 $\nabla^2k(x,y)$ 称为LoG滤波器，表达式为

$\begin{aligned} \nabla^2k(x,y)&=\frac{\partial^2k}{\partial x^2}+\frac{\partial^2k}{\partial y^2}\\ &=\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)\left(\frac1{2\pi\sigma_1\sigma_2}e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)}\right) \end{aligned}$

可以看出 $x,y$ 在函数中是对称的，因此可以先计算

$\begin{aligned} \frac{\partial^2k}{\partial x^2}&=\frac{\partial^2}{\partial x^2}\left(\frac1{2\pi\sigma_1\sigma_2}e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)}\right)\\ &=\frac1{2\pi\sigma_1\sigma_2}\frac\partial{\partial x}\left(-\frac{x}{\pi\sigma_1^2}e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)}\right)\\ &=\frac1{2\pi\sigma_1\sigma_2}\left(\frac{x^2}{\pi\sigma_1^4}-\frac1{\pi\sigma_1^2}\right)e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)} \end{aligned}$

从而可以得出

$\frac{\partial^2k}{\partial y^2}=\frac1{2\pi\sigma_1\sigma_2}\left(\frac{y^2}{\pi\sigma_2^4}-\frac1{\pi\sigma_2^2}\right)e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)}$

则

$\begin{aligned} \nabla^2k(x,y)&=\frac{\partial^2k}{\partial x^2}+\frac{\partial^2k}{\partial y^2}\\ &=\frac1{2\pi^2\sigma_1\sigma_2}\left[\left(\frac{x^2}{\sigma_1^4}-\frac1{\sigma_1^2}\right)e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)}+\left(\frac{y^2}{\sigma_2^4}-\frac1{\sigma_2^2}\right)e^{-\left(\frac{x^2}{2\pi\sigma_1^2}+\frac{y^2}{2\pi\sigma_2^2}\right)}\right] \end{aligned}$

有人问：“博主博主，为什么你的表达式这么复杂，别的博主表达式又很简单呢？”那是因为在OpenCV的高斯滤波中，标准差 $\sigma_1,\sigma_2$ 是可以取不同的值的。特别的，当 $\sigma_1=\sigma_2=\sigma\not=0$ 时，上式变为

$\begin{aligned} \nabla^2k(x,y)&=\frac{\partial^2k}{\partial x^2}+\frac{\partial^2k}{\partial y^2}\\ &=\frac{x^2+y^2-\sigma^2}{2\pi^2\sigma^6}e^{-\frac{x^2+y^2}{2\pi\sigma^2}} \end{aligned}$

关于 $h(x,y)=\nabla^2((f*k)(x,y))=(f*\nabla^2k)(x,y)$ 的证明：首先考虑一维可导函数 $f(x),g(x)$ ，我们有函数的卷积的导数等于一个函数与另一个函数的导数的卷积，即 $\frac{d(f*g)}{dx}(x)=(f*\frac{dg}{dx})(x)$ 。证明如下：

$\begin{aligned} \frac{d(f*g)}{dx}(x)&=\frac{d}{dx}\int_{-\infty}^\infty f(\tau)g(x-\tau)d\tau\\ &=\int_{-\infty}^\infty f(\tau)\frac{d}{dx}g(x-\tau)d\tau\\ &=\left(f*\frac{dg}{dx}\right)(x) \end{aligned}$

同理，考虑二维可导函数 $f(x,y),g(x,y)$ ，我们有函数的卷积的偏导数等于一个函数与另一个函数的偏导数的卷积，即 $\frac{\partial(f*g)}{\partial x}(x,y)=(f*\frac{\partial g}{\partial x})(x,y)$ 。证明如下：

$\begin{aligned} \frac{\partial(f*g)}{\partial x}(x,y)&=\frac{\partial}{\partial x}\int_{-\infty}^\infty\int_{-\infty}^\infty f(\tau,\mu)g(x-\tau,y-\mu)d\tau d\mu\\ &=\int_{-\infty}^\infty\int_{-\infty}^\infty f(\tau,\mu)\frac{\partial}{\partial x}g(x-\tau,y-\mu)d\tau d\mu\\ &=\left(f*\frac{\partial g}{\partial x}\right)(x,y) \end{aligned}$

根据上式，我们很容易证明，二维函数在二阶偏导数的情况。可以用类似上式的方式证明，也可以根据二阶偏导数是一阶偏导数的偏导数来证明：

$\begin{aligned} \frac{\partial^2(f*g)}{\partial^2x}(x,y)&=\frac{\partial}{\partial x}\left(\frac{\partial(f*g)}{\partial x}(x,y)\right)\\ &=\frac{\partial}{\partial x}\left(\left(f*\frac{\partial g}{\partial x}\right)(x,y)\right)\\ &=\left(f*\frac{\partial^2g}{\partial x^2}\right)(x,y) \end{aligned}$

因此我们有

$\begin{aligned} h(x,y)&=\nabla^2((f*k)(x,y))\\ &=\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)(f*k)(x,y)\\ &=\left(f*\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)k\right)(x,y)\\ &=(f*\nabla^2k)(x,y) \end{aligned}$