索伯算子

索伯算子（Sobel operator）是图像处理中的算子之一，有时又称为索伯-费尔德曼算子或索伯滤波器，在影像处理及电脑视觉领域中常被用来做边缘检测。索伯算子最早是由美国计算机科学家艾尔文·索伯及盖瑞·费尔德曼（Gary Feldman）于1968年在史丹佛大学的人工智能实验室（SAIL）所提出，因此为了表扬他们的贡献，而用他们的名字命名。在技术上，它是一离散性差分算子，用来运算图像亮度函数的梯度之近似值。在图像的任何一点使用此算子，索伯算子的运算将会产生对应的梯度向量或是其范数。概念上，索伯算子就是一个小且是整数的滤波器对整张影像在水平及垂直方向上做卷积，因此它所需的运算资源相对较少，另一方面，对于影像中的频率变化较高的地方，它所得的梯度之近似值也比较粗糙。

## 公式

该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即按所示方程式计算：

$$\begin{align} & \left( \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \end{bmatrix} * \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix} \right) [2,2] = \\ & (i \cdot 1)+(h \cdot 2)+(g \cdot 3)+ \\ & (f \cdot 4)+(e \cdot 5)+(d \cdot 6)+ \\ & (c \cdot 7)+(b \cdot 8)+(a \cdot 9) \end{align}$$

即可分别得出横向及纵向的亮度差分近似值。如果以$\mathbf{A}$代表原始图像，$\mathbf{G_x}$及$\mathbf{G_y}$分别代表经横向及纵向边缘检测的图像，其公式如下：

$$\mathbf{G_x} = \begin{bmatrix} +1 & 0 & -1 \\ +2 & 0 & -2 \\ +1 & 0 & -1 \end{bmatrix} * \mathbf{A} \quad \mbox{and} \quad \mathbf{G_y} = \begin{bmatrix} +1 & +2 & +1 \\ 0 & 0 & 0 \\ -1 & -2 & -1 \end{bmatrix} * \mathbf{A}$$

图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大小。

$$\mathbf{G} = \sqrt{ \mathbf{G_x}^2 + \mathbf{G_y}^2 }$$

然后可用以下公式计算梯度方向。

$$\mathbf{\Theta} = \operatorname{arctan}\left({ \mathbf{G_y} \over \mathbf{G_x} }\right)$$

以纵向边缘为例，如果角度$\Theta$等于零，即代表图像的纵向边缘右方较亮；若为 $\pi$，则左方较亮。

## 更正式来说

虽然索伯算子代表着相对没那么精准的影像梯度近似，但其品质也足够在很多应用派上用场。更精确来说，对于每个点，它只用其亮度的值在周围3x3的区域去近似该点的梯度，而且在近似的过程，它也只用整数去代表影像亮度的权重。

## 索伯算子在其他维度的延伸

索伯算子包含着两个分开的运算:

- 利用三角形函数滤波器来达到平滑化 : $h(-1) = 1, h(0) = 2, h(1) = 1$
- 一次微分的中央差值 : $h'(-1) = 1, h'(0) = 0, h'(1) = -1$

索伯算子在不同维度的表达形式，$x,y,z,t \in \left\{0, -1, 1\right\}$:

1D: $h_x'(x) = h'(x);$

2D: $h_x'(x,y) = h'(x)h(y)$

3D: $h_x'(x,y,z) = h'(x)h(y)h(z)$

4D: $h_x'(x,y,z,t) = h'(x)h(y)h(z)h(t)$

下面是一个三维的索伯算子在z轴的实例:

$$h_z'(:,:,-1) = \begin{bmatrix} +1 & +2 & +1 \\ +2 & +4 & +2 \\ +1 & +2 & +1 \end{bmatrix} \quad h_z'(:,:,0) = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad h_z'(:,:,1) = \begin{bmatrix} -1 & -2 & -1 \\ -2 & -4 & -2 \\ -1 & -2 & -1 \end{bmatrix}$$

## 技术细节

索伯算子在实作上可以用很简单的硬件和软件: 只有该点的周围八个点需要被计算，同时在梯度向量的近似上也都只牵涉到整数的计算，而且上述的两个离散滤波器是可分解的:

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

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

因此 Gx 和 Gy 可推导成

$$\mathbf{G}_x = \begin{bmatrix} 1 \\ 2 \\ 1 \end{bmatrix} * \left ( \begin{bmatrix} 1 & 0 & -1 \end{bmatrix} * \mathbf{A} \right ) \quad \mbox{and} \quad \mathbf{G}_y = \begin{bmatrix} \ \ 1 \\ \ \ 0 \\ -1 \end{bmatrix} * \left ( \begin{bmatrix} 1 & 2 & 1 \end{bmatrix} * \mathbf{A} \right )$$

在特别的实作上，可分解的运算可能会有优势因为在每个点的运算上会有少一些的算术运算。

应用卷积核K在像素群P上的虚拟码可以表示成:

N(x,y) = Sum of { K(i,j).P(x-i,y-j)}, for i,j running from -1 to 1.

P代表着像素矩阵，所以N(x,y) 代表的应用卷积核K在像素群P后产生的矩阵

## 实例

经过索伯算子所产生的结果是一张二维的图，代表的每个点的梯度。下图中高梯度的地方会以白色的线表示。为了说明这一点，下方展示了在一张简单的图片上使用索伯算子的结果。

<image>

｜

<image>

<image>

｜

<image>

下方的灰阶影像说明了梯度的方向。红色和黄色代表梯度角度是正的，相对地，蓝色和青色代表梯度角度是负的。在圆的最左边和最右边，垂直的梯度角度为零是因为没有区域的变化，而在圆的最顶部和底部，水平的梯度角度为零是因为没有区域的变化。在圆的最顶部和底部，水平梯度角度变化分别是 −π/2 和π/2是因为没有区域的变化。顶部的负梯度角度代表这条边是从亮到暗的区域界线，而底部的正梯度角度代表这条边是从暗到亮的区域界线。其他显示为黑色的像素是因为附近完全没有任何变化。值得注意的是，因为角度是对应到像素值的函数，因此可能一些微小的变化也可能会导致很大的角度变化。因此，些微的噪声可能会对角度响应变化，而这通常是不乐见的情况。由上述可知，使用角度梯度在影像处理的应用时，需要在影像去噪多加着墨以减少错误。

<image>

｜

<image>

## 替代的运算子

索伯算子在转动上没有完美的对称。因此沙尔想要改进这个特性。它曾提出了更大的5x5核心，不过后来最常被使用的是:

$$\begin{bmatrix} +3 & 0 & -3 \\ +10 & 0 & -10 \\ +3 & 0 & -3 \end{bmatrix} \ \ \ \ \ \ \ \ \ \begin{bmatrix} +3 & +10 & +3 \\ <pre> 0 & 0 & 0 \\ </pre> -3 & -10 & -3 \end{bmatrix}$$

沙尔算子的概念源自于试图在频域上最小化加权的均方角度差，这个最佳化问题需要考量滤波器在数值上的一致性，因此它们是微分核。

另一个相似的改进方法是被法立德和西蒙切利所提出，它们的研究是针对更高次方的微分，相较于沙尔算子，这些滤波器没有考量到数值的一致性。

除此之外，克龙也提出了利用数值方法的最佳化来设计微分滤波器。

2014年，哈斯特也利用任意的三次样条曲线去设计微分滤波器，他的研究指出可以利用长度为7的滤波器作双重滤波在三次样条曲线及三角样条上已取得精准的一次微分和二次微分。

另外一个相似的算子是卡雅利算子，概念也是来自于索伯算子，但它是个具有完美的旋转对称的3x3卷积滤波器。

在那之后的工作，沙尔又针对了光流的预测上提出了更好的滤波器，同一年，他也提出了更好的二次微分滤波器应用在动作预测上。结果显示使用的核心越大，就越有可能去逼近高斯滤波的微分。

## 结果比较

下图是四张不同的梯度算子所预测出的测试图的梯度强度

<image>

｜

<image>

｜

<image>

<image>

｜

<image>

## 伪代码

```
function sobel(A : as two dimensional image array)
	Gx=[-1 0 1; -2 0 2; -1 0 1]
	Gy=[-1 -2 -1; 0 0 0; 1 2 1]

	rows = size(A,1)
	columns = size(A,2)
	mag=zeros(A)

	for i=1:rows-2
		for j=1:columns-2
			S1=sum(sum(Gx.*A(i:i+2,j:j+2)))
			S2=sum(sum(Gy.*A(i:i+2,j:j+2)))

			mag(i+1,j+1)=sqrt(S1.^2+S2.^2)
		end for
	end for

	threshold = 70 %varies for application [0 255]
	output_image = max(mag,threshold)
	output_image(output_image==round(threshold))=0;
	return output_image
end function
```