图像处理中的散度

gljin_cn 2014-07-28

展开全文

1. 什么是散度？
2. 散度在图像处理中有何应用

1. 什么是散度
1.1 散度的定义
散度是作用在向量场上的一个算子。
用三维空间来举例，向量场就是在空间每一点处都对应一个特殊的三维向量的向量函数：
F(x,y,z)=(v_1(x,y,z),v_2(x,y,z),v_3(x,y,z))^T

。比如海洋里，各点在单位时间单位体积中水的流量就是一个三维场，这个场也称为通量。

散度算子定义为：
$div(F)=\frac{\partial v_1}{\partial x}+\frac{\partial v_2}{\partial y}+\frac{\partial v_3}{\partial z}$ 。
它是一个标量函数(场)，也就是说，在定义空间中每一点的散度是一个值。

1.2 散度的物理意义
用水流来解释，散度的物理意义可以叙述为：

如果一点的散度大于0，那么在这一点有一个水龙头不断往外冒水（称为源点）
如果一点的散度小于0，那么在这一点有一个下水道，总有一些水只进不出（称为汇点）
如果一点的散度等于0，那么请放心，在这个点周围的小区域里，单位时间进来多少水就出去多少水。

1.3 数学推导
咱们来看看在一点 A(x,y,z)

的附近到底发生了什么。以这一点为中心，用一个边长分别为 dx,dy,dz

的平行于坐标轴的长方体盒子包围它，来详细分析长方体各表面水向外跑了多少。先看盒子在

方向上的两个面：
第一个面是一个面积为 dydz

的长方形，它的中心坐标是 $(x+\frac{dx}{2},y,z)$ ，这一点的通量是 $F(x+\frac{dx}{2},y,z)$ ，用Taylor展开式可以近似为： $(v_1,v_2,v_3)^T+\frac{dx}{2}(\frac{\partial v_1}{\partial x}, \frac{\partial v_2}{\partial y}, \frac{\partial v_3}{\partial z})^T$ ，又因为这一长方形的外法线方向是 (1,0,0)

，因此这一面在单位时间向外的流量就是二者相乘再乘以面积，由于法线的特殊形式，y、z分量自动消失了：
$V_{x+}=(v_1+\frac{dx}{2}\frac{\partial v_1}{\partial x})dydz$
同理，在x负半轴上的那个面单位时间向外的流量是：
$V_{x-}=(-v_1+\frac{dx}{2}\frac{\partial v_1}{\partial x})dydz$
因此单位时间在x方向上的总的向外的流量是：
$V_x=V_{x+}+V_{x-}=\frac{\partial v_1}{dx}dxdydz$
把三个坐标轴向外的流量加在一起，我们就得到了围绕点

，体积为

的长方体单位时间向外的流量是 $V=\left( \frac{\partial v_1}{\partial x}+\frac{\partial v_2}{\partial y}+\frac{\partial v_3}{\partial z} \right) dxdydz=div(F)dxdydz$ 。

从上面的推导立即可以得出结论：

在一个区域中，单位时间向外的总流量就是把每一个小区域向外的流量加起来(内部相互抵消，最终只有区域边界上的值得以展现)：

$V(D)=\int_{D}div(F)dxdydz$

平均到一个点上，单位时间向外流量的密度就是 $\frac{div(F)dxdydz}{dxdydz}=div(F)$
一个区域无论多复杂，只要不包含源点和汇点，其上散度的积分一定为0

1.4 散度与扩散
假设在空间中有一个浓度场 $\rho (x,y,z)$ ，则在每一点都有一个浓度上升最快的方向，我们称其为梯度 $\nabla \rho =\left( \frac{\partial \rho}{\partial x}, \frac{\partial \rho}{\partial y}, \frac{\partial \rho}{\partial z} \right)$ ，它是一个向量场。浓度差带来的后果就是空间物质会发生运动，从高浓度向低浓度运动，其结果就是浓度中和,趋向平衡。这种运动可以用一个偏微分方程描述：
$\frac{\partial \rho }{\partial t}=-div(K \cdot -\nabla \rho)=div(K \cdot \nabla \rho)$
这个方程被形象地称为扩散方程，来源于物理上的连续性方程。等式右边一定是负散度，因为若一个点散度为正，说明它浓度大，扩散应该减少它的值，然而因为浓度对应的运动场是梯度的负值（高浓度向低浓度流动），因此恰好内外两个负号抵消了，最终右端就出现了貌似不科学的正散度结果，不要被迷惑住。散度算子内部的量

可以是标量，也可以是矩阵，用于调节浓度差与扩散方向之间的关系。

为了更加直观地理解，咱先略去多余因子，这样方程就变成了：
$\frac{\partial \rho }{\partial t}=div(\nabla \rho)=\rho_{xx}+\rho_{yy}+\rho_{zz}=\Delta \rho$
等式右边被称为Laplace算子，一般用一个正三角来简写，你可以用二阶导数来理解它。在一小段时间间隔上，这个方程又可以离散化为：
$\rho^{(t+1)}=\rho^{(t)} + dt \cdot div(\nabla\rho^{(t)})=\rho^{(t)} + dt \cdot \Delta \rho ^{(t)}$
直接含义就是：在每个小时间段内，如果一个点的二阶导数大于0，则把它的浓度增加一些，如果一个点二阶导数小于0，则把它的浓度降低一些。因为二阶导数大于0的点往往是下凹的点，是局部极小值，因此增加它可以让局部浓度变平滑;类似地，二阶导数小于0的点往往是上凸点，是局部极大值，要减少它才能更平滑。

当时间趋向于无穷大时，方程达到稳定，左端为0，那么我们就得到稳定值满足的条件：整个区域上散度为0。也可以理解为最终消灭了所有的源点和汇点，场变得光滑了。

2. 散度在图像去噪中的应用
在图像领域散度算子主要用在去噪中。假设一幅图像为 I(x,y)

，它的梯度算子 $\nabla I=(\frac{\partial I}{\partial x}, \frac{\partial I}{\partial y})^T$ 是一个二维场，那么我们立即可以用散度算子构造一个扩散方程：
$\frac{\partial I}{\partial t}=div(\nabla I)=I_{xx}+I_{yy}$
把这个扩散方程作用于图像就可以去噪了，上面已经解释了它的作用过程是比较图像上的每个点，如果一个点值比周围点低，就增加它，如果比周围点高，就减少它，实质就是平滑图像。但是由于它是各向同性的均匀扩散方程，导致图像上所有细节均匀模糊，去噪效果很糟糕。

Perona和Malik在90年代初发现，由于图像边缘往往处在梯度值较大的点处，如果扩散方程在梯度值较大的区域减速扩散，在梯度值较小的区域加速扩散，则可以在着重去噪的同时保护图像有用细节。他们修改后的扩散方程就是有名的P-M方程：
$\frac{\partial I}{\partial t}=div(g(|\nabla I|)\cdot \nabla I)\approx c_1 I_{\xi \xi }+c_2I_{\eta \eta }$

其中函数g是一个递减函数，保证随图像梯度模值增大函数值递减，起到只在图像平滑区域（小梯度点）猛烈扩散的作用。同时，这个方程还可以变形为在图像局部沿边缘方向 $\xi$ 和跨边缘方向 $\eta$ 上的两个一维扩散之和，好的算法能保证在沿着边缘方向扩散地多，跨边缘扩散地少，也就是保证 c_1>c_2

，起到在去噪的同时保护边缘的作用。散度形式和方向导数形式的扩散方程是随后P-M方程改进的两个主要方向。

基于扩散方程的去噪方法的优点主要有：

结合微分几何和物理方程，比较高大上；
可以控制图像局部区域的扩散特性，对图像的控制力强；
易于推广到三维和更高维以及流形（比如地球表面）上，方程都不用变。

缺点主要有：

速度慢，因为是迭代算法；
扩散会导致边缘发生一定程度地移位；
理论难于往深发展。

最后请欣赏梵高的名画星空：

各向同性扩散方程对其进行均匀平滑的结果（啥都看不清了）：
如果第一幅图是初始浓度，这一幅图就是一段时间后的中和了的浓度。最终，图像成为一张纯色画布，恰好是原始图像的均值。这一过程在数学上严格等价于对图像做方差不断增大的高斯卷积。

修改扩散系数后方向可控平滑的结果（只沿着边缘扩散，保护边缘）：
明显可以看出，在平滑噪声（细小颜色杂质）的同时，大边缘得到了很好地保留。

注：原文关于函数g的叙述有错，感谢xiao huang的细心观察！
参考文献
1. R.P.Feynman et al. 费恩曼物理学讲义（第二卷）.上海科学技术出版社, 2005.
2. 王大凯, 侯榆青,彭进业. 图像处理的偏微分方程方法. 科学出版社, 2008.
3. 王小龙,彭国华. 彩色图像的方向扩散去噪模型研究[J]. 计算机工程与应用, 2013, 49(22): 208-211