Roberts算子、Sobel算子和Laplacian算子的数学推导

温馨提示:本文数学公式推导较多,建议在大屏下(电脑或者平板)观看比较合适)

准备基础知识

欧氏距离(欧几里得度量)

定义

是一个通常采用的距离定义,指在m维空间中两个点之间的真实距离,或者向量的自然长度(即该点到原点的距离)。在二维和三维空间中的欧氏距离就是两点之间的实际距离

计算公式

二维平面上两点间的欧氏距离
$d_{12}=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$

三维空间两点间的欧氏距离
$d_{12}=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}$

两个n维向量间的欧氏距离
$d_{12}=\sum_{k=1}^n(x_{1k}-x_{2k})^2$

城市街区距离(曼哈顿距离、出租车距离)

定义

用来标明两个点在标准坐标系上的绝对轴距总和

举例说明

如下图所示:绿线代表欧式距离红线代表曼哈顿距离黄线和蓝线代表等价的曼哈顿距离
城市街区距离图示

计算公式

二维平面两点间的曼哈顿距离
$d_{12}=|x_1-x_2|+|y_1-y_2|$

两个n维向量间的曼哈顿距离
$d_{12}=\sum_{k=1}^n|x_{1k}-x_{2k}|$

一阶导数

连续函数的情况

对于连续函数,微分表达式为
$f’(x)=\lim\limits_{h\rightarrow0}\frac{f(x+h)-f(x)}{h}$
或者
$f’(x)=\lim\limits_{h\rightarrow0}\frac{f(x+h)-f(x-h)}{2h}$

离散函数的情况

对于离散函数(图像),其导数必须用差分方程来表示,
前向差分
$I_x=\frac{I(x)-I(x-h)}{h}$
后向差分
$I_x=\frac{I(x+h)-I(x)}{h}$
以及中心差分
$I_x=\frac{f(x+h)-f(x-h)}{2h}$

图像梯度

图像I的梯度定义为
$\nabla(I(x,y))=\left[I_x=\frac{\partial I}{\partial x},I_y=\frac{\partial I}{\partial y}\right]$
幅值
$||\nabla(I)||=\sqrt{I_x^2+I_y^2}$
幅值也可用
$||\nabla(I)||\approx|I_x|+|I_y|$
近似

二阶导数

连续函数的情况

对于连续函数(假设其二阶可导)
二阶导数
$f’’(x)=\lim\limits_{h\rightarrow0}\frac{f’(x+h)-f’(x)}{h}$
也就是
$f’’(x)=\lim\limits_{h\rightarrow0}\frac{f(x+h)-2f(x)+f(x-h)}{h^2}$

离散函数的情况

对于离散函数(图像),其差分方程
$I_{xx}=\frac{I(x+h)-2I(x)+I(x-h)}{h^2}$

罗伯茨算子(Roberts)

定义

Roberts算子是一种斜向偏差分的梯度计算方法,梯度的大小代表边缘的强度,梯度的方向与边缘的走向垂直。Roberts操作实际上是求旋转45度两个方向上微分值的和

计算方法

如果我们沿如下图方向角度求其交叉方向的偏导数,则得到Roberts于1963年提出的交叉算子边缘检测方法(这个结论是我根据推导的结果总结得到的,并非预先就知道的)。该方法最大优点是计算量小,速度快。但该方法由于是采用偶数模板,如下图所示,所求的(x,y)点处梯度幅度值,其实是图中交叉点处的值,从而导致在图像(x,y)点所求的梯度幅度值偏移了半个像素(见下图)
罗伯茨算子

$\frac{\partial f}{\partial x}\approx f(x+1,y)-f(x,y)\approx f(x+1,y+1)-f(x,y+1)$
$\frac{\partial f}{\partial y}\approx f(x,y+1)-f(x,y)\approx f(x+1,y+1)-f(x+1,y)$
$||\nabla f(x,y)||=\sqrt{(\frac{\partial f}{\partial x})^2+(\frac{\partial f}{\partial y})^2}\approx |\frac{\partial f}{\partial x}|+|\frac{\partial f}{\partial y}|$
$①\approx \frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}\approx (f(x+1,y)-f(x,y))+(f(x+1,y+1)-f(x+1,y))=f(x+1,y+1)-f(x,y)$
$②\approx \frac{\partial f}{\partial x}-\frac{\partial f}{\partial y}\approx (f(x+1,y)-f(x,y))-(f(x,y+1)-f(x,y))=f(x+1,y)-f(x,y+1)$
$③\approx \frac{\partial f}{\partial y}-\frac{\partial f}{\partial x}\approx (f(x,y+1)-f(x,y))-(f(x+1,y)-f(x,y))=f(x,y+1)-f(x+1,y)$
$④\approx -\frac{\partial f}{\partial x}-\frac{\partial f}{\partial y}\approx -(f(x+1,y)-f(x,y))-(f(x+1,y+1)-f(x+1,y))=f(x,y)-f(x+1,y+1)$

Roberts算子
$G’x=
\left[
\begin{matrix}
{1} & {0} \\
{0} & {-1}
\end{matrix}
\right]
$

$G’y=
\left[
\begin{matrix}
{0} & {1} \\
{-1} & {0}
\end{matrix}
\right]
$

索贝尔算子(Sobel)

定义

Sobel算子考虑了水平、垂直和2个对角共计4个方向对的梯度加权求和,是一个3×3各向异性的梯度算子
注:并不是现行教科书上描述的简单的隔行/列的差分运算,中心像素位置也并未参与运算

数学基础

  • 笛卡尔网格

  • 前向差分

  • 距离反比的4方向对梯度加权

  • 城市街区距离
    注:这些概念已经在上面准备知识部分进行了部分说明

像素N(p)邻域及笛卡尔网格

像素N(p)邻域及笛卡尔网格

计算方法

(1)定义一个给定邻域方向梯度矢量g的幅度为 $|g| = <像素灰度差分>/<相邻像素的距离>$

(2)Sobel采用的像素距离是一种城市街区距离而并非通常的欧式距离。因此,对角方向相邻像素之间的距离值为2

(3)矢量g的方向可以通过中心像素相关邻域的单位矢量给出,这里的邻域是对称出现的,即四个方向对:。沿着4个方向对上求其梯度矢量和,可以给出当前像素平均梯度估计,则有

$G=(z_1-z_9)\cdot\frac{1}{4}\cdot[-1,1]+(z_2-z_8)\cdot\frac{1}{2}\cdot[0,1]+(z_3-z_7)\cdot\frac{1}{4}\cdot[1,1]+(z_6-z_4)\cdot\frac{1}{2}\cdot[1,0]$
$=[(z_3-z_7-z_1+z_9)\cdot\frac{1}{4}+(z_6-z_4)\cdot\frac{1}{2},(z_3-z_7+ z_1-z_9)\cdot\frac{1}{4} + (z_2-z_8)\cdot\frac{1}{2}]$

式中, 4个单位向量 [1, 1],[-1, 1],[0, 1], [1, 0] 控制差分的方向,系数1/4, 1/2为距离反比权重
(4)理论上应该对G除以4得到平均梯度估计除法会丢失低阶的重要字节, 更方便的是把向量乘于4,而不是除于4,以保留低阶字节。因此,计算出的估计值比平均梯度在数值上扩大了16倍

$G’=4\cdot G$

$=[z_3-z_7-z_1+z_9+2\cdot (z_6-z_4),z_3-z_7+z_1-z9+2\cdot (z_2-z_8)]$

$=[z_3+2\cdot z_6+z_9-z_1-2\cdot z_4-z_7,z_1+2\cdot z_2+z_3-z_7-2\cdot z_8-z_9]$

(5)按x-y方向,可分别写成:

$G’x = (z_3 + 2\cdot z_6 + z_9) - (z_1 + 2\cdot z_4 + z_7)$

$G’y = (z_1 + 2\cdot z_2 + z_3) - (z_7 + 2\cdot z_8 + z_9)$

即Sobel的横向和纵向的滤波矩阵
$G’x=
\left[
\begin{matrix}
{-1} & {0} & {1} \\
{-2} & {0} &{2} \\
{-1} & {0} & {1}
\end{matrix}
\right]
$

$G’y=
\left[
\begin{matrix}
{1} & {2} & {1} \\
{0} & {0} &{0} \\
{-1} & {-2} & {-1}
\end{matrix}
\right]
$

拉普拉斯算子(Laplacian)

对于连续函数的推导

拉普拉斯算子是n维欧式空间的一个二阶微分算子,其定义为两个梯度向量算子的内积(也就是梯度的散度)
$\Delta=\nabla\cdot\nabla$
$=\left(\frac{\partial }{\partial x}\vec{i}+\frac{\partial }{\partial y}\vec{j}\right)\cdot\left(\frac{\partial }{\partial x}\vec{i}+\frac{\partial }{\partial y}\vec{j}\right)$
$=\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}$
其在二维空间上的表达式为
$\Delta f=\frac{\partial f^2 }{\partial x^2}+\frac{\partial f^2 }{\partial y^2}$

对于离散函数的推导

注:对于离散函数而言,导数退化为差分,因为相邻点像素距离均为1,故下面分母上的1均被省略

一维离散的情况

对于一维离散来说,
一阶差分
$\frac{\partial f}{\partial x}=f(x+1)-f(x)$
二阶差分
$\frac{\partial^2 f}{\partial x^2}=f(x+1)-2f(x)+ f(x-1)$

因此,1维拉普拉斯运算可以通过1维卷积核[1,-2,1]实现

二维离散的情况

对于二维离散来说,f(x,y)对x右侧的一阶差分
$\frac{\partial f}{\partial x} = f(x+1,y)-f(x,y)$
f(x,y)对x左侧的一阶差分
$\frac{\partial f}{\partial x} = f(x,y)-f(x-1,y)$
f(x,y)对x的二阶差分(右侧一阶减左侧一阶,分母仍然是1略去)
$\frac{\partial^2 f}{\partial x^2}=(f(x+1,y)-f(x,y))-(f(x,y)-f(x-1,y))$
$=f(x+1,y)+f(x-1,y)-2f(x,y)$
同理对f(x,y)对y的二阶差分
$\frac{\partial^2 f}{\partial y^2}= f(x,y+1)-f(x,y)-(f(x,y)-f(x,y-1))$
$=f(x,y+1)+f(x,y-1)-2f(x,y)$
将上述两式子相加得(注:拉普拉斯算子是2个维上二阶差分的和
$\bigtriangledown^2 f= f(x+1,y)+ f(x-1,y)-2f(x,y)+f(x,y+1)+f(x,y-1)-2f(x,y)$
$=(f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1))-4f(x,y)$

因此,2维拉普拉斯运算可以通过以下2维卷积核实现
$
\left[
\begin{matrix}
{0} & {1} & {0} \\
{1} & {-4} &{1} \\
{0} & {1} & {0}
\end{matrix}
\right]
$