从边缘检测聊起

从边缘检测聊起

边缘检测的诞生来自两个方面:

  • 一方面,生物学家对人体视网膜的研究,发现人类对边缘线条的变化非常敏感。

  • 另一方面,基于傅里叶变换的一维数字信号处理的快速发展。我们也从一维的时间频率推广到二维的空间频率。根据经验,一维时间信号中的噪声信号或者突变信号往往对应着频域的高频成分。类似地,二维图片的边缘信息对应着空间频域的高频成分。
    如下图所示,在频域的低通处理会去掉原始图片的边缘,换句话说,使图片变得更柔和,反之,高通滤波则会只保留图片的边缘细节。
    具体可以看我的博客: 二维傅里叶变换

由此引出本文的第一个问题

1. 如何获取边缘特征 ?

回答就是计算 梯度 (gradient) 或者说一阶梯度,所谓边缘,就是像素值剧烈变化的地方,而像素的变换,用数学术语来讲,就是一阶梯度值啦。如下图所示,黑白分隔的边缘对应的就是一阶梯度的极大或极小值。

从频域分析,边缘对应着高频成分;从梯度分析,边缘对应着一阶梯度的极值点 (局部极大值,或极小值),或者说二阶梯度的过零点

1.1 梯度 – 差分卷积核

连续梯度计算如下: \[
\cfrac {df} {dx} = \lim_{\Delta x \to 0} \cfrac {f(x) – f(x – \Delta x)} {\Delta x} = f'(x)
\]
对于一维离散信号,令上面的 \(\Delta x = 1\), 可以得到离散梯度: \[
\begin {align}
\text {Backward} \quad \quad \quad &\cfrac {df} {dx} = f(x) – f(x-1) \\
\text {Forward} \quad \quad \quad &\cfrac {df} {dx} = f(x) – f(x+1) \\
\text {Central} \quad \quad \quad &\cfrac {df} {dx} = f(x+1) – f(x-1)
\end {align}
\]
如果我们要计算所有离散点的梯度 \(\cfrac {df} {dx}\) ,那么等价于 原始信号和差分卷积核 (differentiation kernel) 进行卷积操作 (注意需要翻转),如下是一维差分卷积核 \(h\) \[
\begin {align}
\text {Backward} \quad \quad \quad &h_{back} = \begin {bmatrix} 0 & 1 & -1 \end{bmatrix} \\
\text {Forward} \quad \quad \quad &h_{for} =\begin {bmatrix} -1 & 1 & 0 \end{bmatrix} \\
\text {Central} \quad \quad \quad &h_{cent} =\begin {bmatrix} 1 & 0 & -1 \end{bmatrix}
\end {align}
\]
卷积结果 \(y = x \ast h\) 就是每个点对应的一阶离散梯度

那么如何求二维信号的梯度呢 ?

其实我在一篇博客中讲过矩阵的偏导计算,有兴趣的可以看看。

这里直接给出结论,对于二维函数 \(f(x, y)\) ,我们有: \[
\begin {align}
\text {Gradient vector} \quad \quad \quad &\nabla f(x, y) = \begin {bmatrix} \frac {\partial f(x,y) } {\partial x} \\ \frac {\partial f(x,y) } {\partial y}\end {bmatrix} = \begin {bmatrix} f_x \\ f_y\end {bmatrix} \\
\text {Gradient magnitude} \quad \quad \quad &\vert \nabla f(x,y) \vert = \sqrt {f_x^2 + f_y^2} \\
\text {Gradient direction} \quad \quad \quad &\theta = \tan^{-1} (f_y/f_x)
\end {align}
\]

参考 一维信号的卷积核,我们可以得到如下的二维差分卷积核 \(K_x\)\(K_y\)
\[
K_x = \cfrac {1}{3}
\begin {bmatrix}
1 & 0 & -1 \\
1 & 0 & -1 \\
1 & 0 & -1
\end{bmatrix} \quad \quad \quad
K_y = \cfrac {1}{3}
\begin {bmatrix}
1 & 1 & 1 \\
0 & 0 & 0 \\
-1 & -1 & -1
\end{bmatrix}
\]
\(K_x\) 计算水平梯度,检测垂直边缘 – \(K_y\) 计算竖直梯度,检测水平边缘

这里有点绕口,但很好理解。因为梯度方向是下降最快的方向,那么垂直于梯度方向自然就是边缘了。因为边缘线上的值基本是相等的,但边缘两侧的值变化很大,梯度方向自然垂直于梯度。

还是请注意: 卷积不等于相关,所以不要认为上面的值位置颠倒了。

值得一提的是,一个更有效的梯度卷积核是 Sobel kernel,如下所示: \[
K_x = \cfrac {1}{3}
\begin {bmatrix}
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1
\end{bmatrix} \quad \quad \quad
K_y = \cfrac {1}{3}
\begin {bmatrix}
1 & 2 & 1 \\
0 & 0 & 0 \\
-1 & -2 & -1
\end{bmatrix}
\]
它的精髓就在于引入了 Gaussian kernel 用于平滑预处理,消除噪声

1.2 什么是 Gaussian kernel ?

高斯函数真的是非常重要的函数,具体多牛逼大家去看维基百科吧。

在图像处理领域,Gaussian kernel 广泛用于图片模糊 (blur),还有根据 Gaussian kernel 引出的 Laplacian of Gaussian kernel,用于小波变换的 Gabor kernel,还有基于尺度不变性的差分高斯金字塔 (DoG, Difference of Gaussian Pyramid) 模型等等,这些以后再说,其中尺度不变性又是另一个很长的故事,它还和热扩散方程扯上了关系 :joy: 。

不得不提的是,高斯函数的傅里叶变换仍然是高斯函数 !!!

由于高斯函数是连续的,离散卷积核要尽量逼近高斯函数的形状。

下面简单列出了一维和二维的大小等于 3 的高斯卷积核。 \[
\begin {align}
\text {1D Gaussian kernel :} &\begin{bmatrix} 1 & 2 & 1 \end {bmatrix} \\
\text {2D Gaussian kernel :}
&\cfrac {1}{16}\begin{bmatrix}
1 & 2 & 1 \\
2 & 4 & 2 \\
1 & 2 & 1
\end {bmatrix}
\end {align}
\]
OpenCV, skimage 等图像处理库都有根据 kernel size 生成 2D Gaussian kernel 的函数。

为什么 Gaussian kernel 可以用于图片模糊平滑处理 ?

答: Gaussian kernel 是低通滤波器

根据空间域的卷积运算对应频域的乘积运算,可以得到原始图片和高斯核卷积操作,其实在频域上等于原始图片的频谱和高斯核的频谱相乘。而高斯核的频谱仍然是高斯函数,且中心点是频率为 0 的原点。因此两者相乘,只会保留低频成分,而过滤掉高频成分,即高斯核是低通滤波器。

连续的二维高斯核 \[
G_\sigma (x, y) = \cfrac {1} {2\pi \sigma^2} \exp(- \cfrac {x^2 + y^2}{2 \sigma^2})
\]
而离散的 \((2k+1) \times (2k+1)\) 的高斯核,它的每个 entry 的值如下:

\[
h_{ij} = \cfrac {1} {2\pi \sigma^2} \exp \Big(-\cfrac {(i-k)^2+(j-k)^2}{2 \sigma^2} \Big ), 0\leq i,j < 2k+1
\]
Note: 高斯核尺寸越大,标准差 \(\sigma\) 越大,则低通效果越强,图片越模糊

1.3 噪声消除

回到 Sobel kernel ,为什么说 Sobel kernel 的精髓在于引入 Gaussian kernel ?

\(K_x\) 水平梯度算子为例,如下所示: \[
K_x = \cfrac {1}{3}
\begin {bmatrix}
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1
\end{bmatrix} = \cfrac {1}{3} \begin {bmatrix}
1 \\ 2 \\1 \end{bmatrix} \begin {bmatrix}
1 & 0 & -1 \end{bmatrix}
\]
可以看到,\(K_x\) 可以分解为 一维高斯核一维差分核 的外积。

前面说过,Sobel kernel 的目的是为了提取边缘特征 (高频),而 Gaussian kernel 的目的是模糊边缘 (去高频)。这两者不是互相矛盾吗 ?

这里就不得不提到另一个图像处理常见的问题 – 噪声 (noise)

实际的图像处理中,我们通常都不会直接对原始图片进行边缘检测。一个重要原因就是噪声的干扰,而噪声通常是高频信号。

请思考为什么噪声通常是高频信号 ?

答: 从直观上看,噪声像素点看起来和周围的像素值明显不同,不然我们为什么称之为噪声呢。自然地,噪声的梯度值就会很大,于是有限差分滤波器 (Finite difference filter) 对这些噪声点的响应值就会很大,结果通常不是很理想。

常见的噪声有 3 种:

  • 椒盐噪声 (salt and pepper noise): 随机的黑白像素点
  • 脉冲噪声 (impulse noise): 随机的纯白像素点
  • 高斯噪声 (gaussian noise): 白噪声的亮度分布为高斯分布

所以我们要引入 smooth kernel,比如 mean kernel, median kernel,还有 gaussian kernel,它们都是低通滤波器

一个具体的例子,下图中,我们直接对一维信号和差分核卷积,由于噪声的干扰,我们无法看到一阶梯度的极大和极小点,都湮没在噪声信号中了。

现在我们将原始信号和高斯核卷积,效果马上就出来了

1.4 sobel kernel 的具体由来

卷积操作有一个重要性质: \[
\cfrac {d} {dt} [f(t) \ast h(t)] = f(t) \ast \cfrac {d} {dt}h(t)
\]
也就是说,卷积之后求导,等价于先对其中某一个信号求导后再进行卷积操作

证明: \[
\cfrac {d} {dt} [f(t) \ast h(t)] = \cfrac {d} {dt} \int f(\tau) h(t – \tau) d\tau = \int f(\tau) \cfrac {d} {dt}h(t – \tau) d\tau = f(t) \ast \cfrac {d} {dt}h(t)
\]

好的,我们现在知道了,原先我们先用高斯核平滑滤波,然后用一阶差分核得到 \(x, y\) 方向的一阶梯度,现在可以这么做,先对高斯核用一阶差分核得到高斯核的一阶梯度,然后在用它和我们要处理的图像直接卷积,结果是先相同的。

高斯核的一阶梯度是什么核呢 ?

答: Sobel kernel.

如下图所示,我们将离散化的二维高斯核和 一阶差分核 \(\begin {bmatrix} 1 & 0 & -1 \end {bmatrix}\) 卷积得到高斯核 \(x\) 方向的梯度,如下图所示

同理可以得到 \(y\) 方向的梯度,\(x,y\) 方向的高斯一阶梯度如下图所示:

如果我们将其用 \(3 \times 3\) 的 卷积核表征,就是我们之前提到的 sobel kernle 啊。

这和之前提到的 sobel kernle 是 高斯核和差分核外积计算得到 的结论是相同的。 \[
K_x = \cfrac {1}{3}
\begin {bmatrix}
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1
\end{bmatrix} = \cfrac {1}{3} \begin {bmatrix}
1 \\ 2 \\1 \end{bmatrix} \begin {bmatrix}
1 & 0 & -1 \end{bmatrix}
\]

1.5 Laplacian operator

之前我们采用的是一阶梯度卷积核 用于边缘检测,在处理之后的图片中,边缘的亮度在局部是最大的,也就是说边缘对应的梯度值是局部极大或极小值。显然,一阶梯度的极大和极小,还对应着二阶梯度的过零点 (zero-crossing),于是自然引出二阶梯度算子。

用来求解二阶梯度的卷积核就是 Laplacian kernel

一阶梯度有 \(x, y\) 两个方向的卷积核,而拉普拉斯算子(Laplacian operator) 是标量运算符,定义为两个梯度向量的内积,简单点说,就是计算各个方向的二阶梯度之和: \[
\nabla^2 = \begin {bmatrix} \cfrac {\partial} {\partial x_1} & \cdots & \cfrac {\partial} {\partial x_N} \end {bmatrix} \begin {bmatrix} \cfrac {\partial} {\partial x_1} \\ \vdots \\ \cfrac {\partial} {\partial x_N} \end {bmatrix} = \sum_{i=1}^N \cfrac {\partial^2} {\partial x^2_i}
\]
Note:

上面表达式千万别理解为一阶梯度的平方和,而应该是二阶梯度之和

\(N=2\) 的二维空间,我们有 \[
\nabla^2f(x, y) = \cfrac {\partial^2 f} {\partial x^2} + \cfrac {\partial^2 f} {\partial y^2}
\]
Laplacian operator 的特点是: 各向同性,对任何走向的边界和线条进行锐化,无方向性 (undirected),也就是说对图片旋转不敏感

前面提到一阶一维离散梯度: \(\nabla f[n] = f[n+1] – f[n]\)

于是二阶离散梯度: \[
\begin {align}
\nabla^2 f[n] = \nabla (\nabla f[n]) &= \nabla f[n] – \nabla f[n-1] \\
&= (f[n+1] – f[n]) – (f[n] – f[n-1]) \\
&= f[n+1] – 2f[n] + f[n-1]
\end {align}
\]
二阶二维离散梯度等于各自维度二阶梯度之和: \[
\begin {align}
\nabla^2f[m, n] &= \nabla^2_m f[m, n] + \nabla^2_n f[m, n] \\
&= (f[m+1, n] – f[m, n]) – (f[m, n] – f[m-1, n]) \\
&\qquad + (f[m, n+1] – f[m, n]) – (f[m, n] – f[m, n+1]) \\
&= f[m+1,n] + f[m-1, n] + f[m, n+1] + f[m, n-1] – 4f[m, n]
\end {align}
\]
\[
\text {Laplacian kernel} = \begin{bmatrix}
0 & 1 & 0 \\
1 & -4 & 1 \\
0 & 1 & 0
\end{bmatrix}
\]
如果考虑每个点周围 8 个点,可以得到另一个算子: \[
\text {Laplacian kernel} = \begin{bmatrix}
1 & 1 & 1 \\
1 & -8 & 1 \\
1 & 1 & 1
\end{bmatrix}
\]
Note: 有的博客或者论文中 Laplacian kernel 的符号和上面正好完全相反,这也是合理的。

1.6 Laplacian of Gaussian kernel (LoG)

既然 高斯核可以求一阶梯度,自然也是可以求 laplacien operator 的,我们将其称为 Laplacian of Gaussian kernel (LoG) \[
\nabla^2 [G_\sigma (x, y) \ast f(x, y)] = [\nabla^2 G_\sigma (x, y) ] \ast f(x,y) = LoG \ast f(x,y)
\]
其中高斯核 \(G_\sigma (x, y) = \cfrac {1} {2\pi \sigma^2} \exp(- \cfrac {x^2 + y^2}{2 \sigma^2})\)

连续高斯函数求 拉普拉斯算子 (二阶梯度之和),得到: \[
LoG(x, y)= \cfrac {\partial^2 } {\partial x^2} G_\sigma (x, y) + \cfrac {\partial^2} {\partial y^2} G_\sigma (x, y) = (\cfrac {x^2+y^2} {\sigma^4} – \cfrac {2} {\sigma^2}) G_\sigma (x, y)
\]

\(5 \times 5\) 的 LoG kernel 来表征: \[
\begin{bmatrix}
0 & 0 & 1 & 0 & 0 \\
0 & 1 & 2 & 1 & 0 \\
1 & 2 & -16 & 2 & 1 \\
0 & 1 & 2 & 1 & 0 \\
0 & 0 & 1 & 0 & 0
\end{bmatrix}
\]
用 LoG kernel 提取边缘的操作如下:

  1. 将 Log kernel 和原始图片卷积
  2. 提取生成图片中值为 0 的点 (zero-crossings)
  3. 通过阈值筛选过零点,只保留那些比较强的点 (large difference between the positive maximum and the negative minimum),目的是进一步消除噪声带来的弱过零点。

2. 边缘检测的巅峰 – Canny 边缘检测

Sobel kernel 成也 高斯核,败也 高斯核,为什么这么说 ?

固然引入高斯核可以减少噪声,但它对边缘的平滑也是影响非常大的,由此带来的问题就是: 边缘过粗

Canny 提出了一个好的边缘检测器的标准:

  • Good detection: the optimal detector must minimize the probability of false positives (detecting spurious edges caused by noise), as well as that of false negatives (missing real edges)
  • Good localization: the edges detected must be as close as possible to the true edges
  • Single response: the detector must return one point only for each true edge point; that is, minimize the number of local maxima around the true edge.

2.1 Canny 边缘检测

步骤如下:

  1. Gaussian kernel 进行图片平滑

  2. Sobel kernel 或者其他一阶梯度卷积核 计算

  3. 计算 \(G = \sqrt {G_x^2 + G_y^2}\) 以及 \(\theta = \arctan (\frac {G_y}{G_x})\)

  4. non-maxima suppression

    如图所示,将周围 8 个点分为 4 组: 左右,上下,左上和右下,以及右上和坐下
    根据上面获得的 \(\theta\) 值,将中心点的值和其中一组进行比较,如果值比其他两个值都大,则保留;否则,设为 0。
    这里有一个很大的坑,千万千万要注意:
    不管是矩阵表示还是图片存储,我们都规定左上角为坐标原点,在图片中,横轴往右是宽度方向,竖轴往下为高度方向,一般用 \((x, y)​\) 表示一个像素点距离原点的宽度和长度,这也符合笛卡尔坐标轴的坐标表述
    但是,在矩阵表述中,我们一般用 \([row, column]\) 表示一点的所在的行和列,所以两者的顺序刚好反了
    这还不是最坑的,最坑的是,由于竖轴正方向是往下的,所以上面这幅图应该沿 \(x\) 轴翻转 \(180^{\circ}\), 即 顺时针表示 \(0^{\circ} \to 360^{\circ}\) ,而不是上面的逆时针表示,所以千万小心。

  5. double-threshold
    上面更新得到的 梯度矩阵 \(G\) , 将其和两个阈值比较,得到两个矩阵 Strong edges 和 weak edges:

    • 如果 \(G\) 的点大于 high threshold, 那么它在 strong edges 矩阵对应的位置设置为 1, 表明该点是边缘的可能性非常大
    • 如果 \(G\) 的点小于 high threshold, 但大于 low threshold, 那么它在 weak edges 矩阵对应的位置设置为 1,表明该点可能是边缘,也可能是噪声。

    这里有一个经验:
    如果一些长的边缘被切割为短线段,说明 weak edges 被删除过多,建议减小 low threshold;反之,如果 杂散边缘 (spurious edges) 过多,则建议删除增大 high threshold。

  6. link-edge
    Canny 提出上面的 weak edges 矩阵的非零点如果和 strong edges 中的某个非零点相邻,那么认为这个点也是边缘点

  7. 最终得到的 edges 就是我们需要的边缘了

发表评论

电子邮件地址不会被公开。 必填项已用*标注

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d 博主赞过: