图像分析：边缘检测中的梯度算子

0 前言

边缘检测是基于边界的图像分割方法的第一步，边缘就是两个不同的相邻区域之间灰度值不连续或者突变的地方。因此，检测边缘就是，检测灰度明显变化的地方。而边缘位置处灰度的明显变化是可以借助计算灰度的微分来检测的。一般使用一阶微分和二阶微分检测边缘，在边缘位置，一阶微分的幅度值会有局部极值，二阶微分的幅度值会出现过零点。本文主要介绍边缘检测中的一阶微分算子----梯度算子，包括Roberts、Prewitt和Sobel三种算子。

1 基本概念

1.1 图像显示方式（图像坐标系）

想要计算梯度图，就要设计模板卷积，首先要搞明白图像在计算时的坐标系，很多博文对应的模板和坐标系都不匹配，我们在后面的卷积操作中主要使用计算坐标系。

（1）计算坐标系

左图Cameraman所用的坐标系统，常用在图像计算中。其坐标原点在左上角，x轴是水平的，并且向右延伸；y是垂直的，并且向下延伸。 $f(x,y)$ 既可以代表这幅图像，也可以表示 $(x,y)$ 坐标处像素的值。

（2）显示坐标系

右图Lena的坐标系统，常用在屏幕显示中，因为屏幕扫描是从左向右，从上向下进行的，原点在图像的左上角，纵轴标记图像的行，横轴标记图像的列。 $I(r,c)$ 既可以代表这个图像，也可以代表 $(r,c)$ 行列交点处的图像值。

图像的坐标系

1.2 梯度

首先我们要知道的是梯度是一个向量，向量的话有方向和大小，梯度方向指向函数变化最快的方向，大小就是它的模，也是最大的变化率。对于二元函数 $z=f(x,y)$ ，它在点 $(x,y)$ 的梯度就是 $grad f(x,y)$ , 或者 $\nabla f(x,y)$ ，就是：
$grad f (x,y)=\nabla f(x,y)=(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})=f_{x}(x,y)\vec{i}+f_{y}(x,y)\vec{j}$

其中， $\nabla =\frac{\partial}{\partial x} \vec{i}+\frac{\partial}{\partial y} \vec{j}$ ，这个梯度向量的幅度和方向角为
$mag(\nabla f)=\|\nabla f_{(2)} \|=\sqrt {G^{2}_{x}+G^2_{y}}$ $\phi (x, y)=arctan(\frac{G_{y}}{ G_{x}})$

1.3 图像的梯度

下图展示了一个灰度图的数学化表达，像素点 $(x,y)$ 的灰度值是 $f(x,y)$ ,它有八个邻域。

$f(x-1,y+1)$	$f(x,y+1)$	$f(x+1,y+1)$
$f(x-1,y)$	$f(x,y)$	$f(x+1,y)$
$f(x-1,y-1)$	$f(x,y-1)$	$f(x+1,y-1)$

图像在点 $(x,y)$ 的梯度为
$G(x,y) =(G_{x},G_{y})=( f_{x}(x,y),f_{y}(x,y))$
其中 $f_{x}(x,y) = f(x+1, y)-f(x,y)$
$f_{y}(x,y)=f(x,y+1)-f(x,y)$
$f_{x}(x,y)$ 即 $g_{x}$ 对应图像的水平方向， $f_{y}(x,y)$ 即 $g_{x}$ 对应水图像的竖直方向。

1.4 模板卷积

要理解梯度图的生成，就要先了解模板卷积的过程。
模板卷积是模板运算的一种方式，其步骤如下：
（1）将模板在输入图像中漫游，并将模板中心与图像中某个像素位置重合；
（2）将模板上各个系数与模板下各对应像素的灰度相乘；
（3）将所有乘积相加（为保持灰度范围，常将结果再除以模板系数之和，后面梯度算子模板和为0的话就不需要除了）；
（4）将上述运算结果（模板的响应输出）赋给输出图像中对应模板中心位置的像素。

模板卷积运算过程

1.5 梯度图

（1）梯度图的生成和模板卷积有什么不同呢？

其实梯度图生成前面和模板卷积相同，不同的是要生成梯度图，还需要在模板卷积完成后计算在点 $(x,y)$ 梯度的幅值，将幅值作为像素值，这样才算完。 $\color{red}{注意：梯度图上每个像素点的灰度值就是梯度向量的幅度！！！}$ 。

（2）生成梯度图的模板

下图是生成梯度图用到的的水平模板和竖直模板：

最简单的水平梯度模板

最简单的竖直梯度模板

例如，如果只想生成水平方向的梯度图，那么只计算水平方向的梯度 $G_{x}$ ，则梯度图上对应点 $(x,y)$ 的灰度值就是 $g(x,y)=\sqrt{G^{2}_{x}}=|G_{x}|$

如果使用上图的水平方向模板那么点 $(x,y)$ 处的灰度值就是 $g(x,y)=|G_(x)|=|f(x+1,y)-f(x,y)|$
如果使用上图的竖直方向模板那么点 $(x,y)$ 处的灰度值就是 $g(x,y)=|G_(y)|=|f(x,y+1)-f(x,y)|$

（3）实际计算中梯度图的生成

一般是水平方向的 $G_{x}$ 和竖直方向的 $G_{y}$ 各用一个模板，然后结合，那么得到梯度图在点 $(x,y)$ 的灰度值就是
$g(x,y)=\sqrt {G^{2}_{x}+G^2_{y}}$
它就是我们上面说到的梯度的幅值，是以计算以2为范数，对应欧式距离，由于涉及平方和开方运算，计算量比较大。（怎么简化计算呢？？换一种近似计算方式吧！！！）
在真实的梯度图输出计算中，采用以1为范数（对应城区距离）的简单计算方式，即
$g(x,y) = |G_{x}|+|G_{y}|$
另一种简单的方式是以 $\infty$ 为范数（对应棋盘距离），即
$g(x,y) =max( |G_{x}|,|G_{y}|)$

2 梯度算子

2.1 什么是梯度算子呢？

梯度算子是一阶导数算子，像拉普拉斯算子就不是梯度算子，因为它是一个二阶导数算子。
梯度算子其实就是 $G_{x}$ 和 $G_{y}$ 对应模板的组合，也就是说梯度算子包含两个模板(甚至是两个以上，如果还有对角线模板)！！

2.2 常用的梯度算子有哪些呢？

首先了解下梯度算子的设计，一般是水平方向和竖直方向，水平方向模板转置再对折就是竖直方向。

(1) Roberts交叉算子

Roberts交叉算子

其本质是一个对角线方向的梯度算子，对应的水平方向和竖直方向的梯度分别为 $G_{x} = f(x+1,y+1)-f(x,y)$ $G_{y} = f(x,y+1)-f(x+1,y)$
输出梯度图在 $(x,y)$ 的灰度值为 $g(x,y)=\sqrt{G^{2}_{x}+G^{2}_{y}}$
优点：边缘定位较准
缺点：（1）没有描述水平和竖直方向的灰度变化，只关注了对角线方向，容易造成遗漏。（2）鲁棒性差。由于 $(x,y)$ 点本身参加了梯度计算，不能有效的抑制噪声的干扰。
适用于边缘明显且噪声较少的图像。

(2) Prewitt算子

Prewitt算子

Prewitt算子是典型的 $3\times3$ 模板，其模板中心对应要求梯度的原图像坐标 $(x,y)$ ， $(x,y)$ 对应的8-邻域的像素灰度值如下表所示：

$f(x-1,y+1)$	$f(x,y+1)$	$f(x+1,y+1)$
$f(x-1,y)$	$f(x,y)$	$f(x+1,y)$
$f(x-1,y-1)$	$f(x,y-1)$	$f(x+1,y-1)$

通过Prewitt算子的水平模板 $M_{x}$ 卷积后，对应的水平方向梯度为
$G_{x} = f(x+1,y+1)-f(x-1,y+1)+f(x+1,y)-f(x-1,y)+f(x+1,y-1)-f(x-1,y-1)$
通过Prewitt算子的竖直模板 $M_{y}$ 卷积后，对应的竖直方向梯度为
$G_{y} = f(x-1,y+1)-f(x-1,y-1)+f(x,y+1)-f(x,y-1)+f(x+1,y+1)-f(x+1,y-1)$
输出梯度图在 $(x,y)$ 的灰度值为
$g(x,y)=\sqrt{G^{2}_{x}+G^{2}_{y}}$
Prewitt算子引入了类似局部平均的运算，对噪声具有平滑作用，较Roberts算子更能抑制噪声。

(3) Sobel算子

Sobel算子

Sobel算子其实就是是增加了权重系数的Prewitt算子，其模板中心对应要求梯度的原图像坐标，对应的8-邻域的像素灰度值如下表所示：

$f(x-1,y+1)$	$f(x,y+1)$	$f(x+1,y+1)$
$f(x-1,y)$	$f(x,y)$	$f(x+1,y)$
$f(x-1,y-1)$	$f(x,y-1)$	$f(x+1,y-1)$

通过Prewitt算子的水平模板 $M_{x}$ 卷积后，对应的水平方向梯度为
$G_{x} = f(x+1,y+1)-f(x-1,y+1)+2f(x+1,y)-2f(x-1,y)+f(x+1,y-1)-f(x-1,y-1)$
通过Prewitt算子的竖直模板 $M_{y}$ 卷积后，对应的竖直方向梯度为
$G_{y} = f(x-1,y+1)-f(x-1,y-1)+2f(x,y+1)-2f(x,y-1)+f(x+1,y+1)-f(x+1,y-1)$
输出梯度图在 $(x,y)$ 的灰度值为
$g(x,y)=\sqrt{G^{2}_{x}+G^{2}_{y}}$
Sobel算子引入了类似局部加权平均的运算，对边缘的定位比要比Prewitt算子好。

3 代码实现

Python调用OpenCV接口实现Sobel算子边缘检测

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt

img = [cv.imread]('dave.jpg',0)
laplacian = [cv.Laplacian](img,cv.CV_64F)
sobelx = [cv.Sobel](img,cv.CV_64F,1,0,ksize=5)
sobely = [cv.Sobel](img,cv.CV_64F,0,1,ksize=5)

plt.subplot(2,2,1),plt.imshow(img,cmap = 'gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(2,2,2),plt.imshow(laplacian,cmap = 'gray')
plt.title('Laplacian'), plt.xticks([]), plt.yticks([])
plt.subplot(2,2,3),plt.imshow(sobelx,cmap = 'gray')
plt.title('Sobel X'), plt.xticks([]), plt.yticks([])
plt.subplot(2,2,4),plt.imshow(sobely,cmap = 'gray')
plt.title('Sobel Y'), plt.xticks([]), plt.yticks([])
plt.show()