各向异性滤波记录

通常用在图像去噪上，文章许多文字来自别的文章，以下是参考链接：

1. https://zhuanlan.zhihu.com/p/443182644
2. https://www.zhihu.com/question/430771515/answer/2738003246
3. https://www.zhihu.com/question/430771515/answer/1597682592

基本思想

传统滤波时，各个方向都一视同仁，边缘和平坦区域都做平滑，去噪的同时容易丢失边缘等有意义的高频，如高斯滤波、均值滤波等。

因此考虑在做滤波时，各个方向不是一视同仁的，边缘等高频跟平坦区域会区别对待。其实双边滤波也是类似的思路。

公式解释

Perona P, Malik J 于 1990 年提出了 PM 扩散方程，其基本公式为

\[ \frac{\partial u}{\partial t} = \text{div}(c(u)\nabla u) \]

其中 $u$ 是图像，$t$ 是第几次迭代，$c(u)$ 是扩散的传导函数，$\nabla u$ 是梯度，$\text{div}$ 是散度，有些会在下面详细解释。

迭代

所以每次迭代，计算 $u(t+1) = u(t) + \lambda * \partial u / \partial t$ 进行更新就行。

这里的 $\lambda$ 取值越大越平滑。

散度计算

散度的这个概念是高数后面开始涉及的，就是高斯-斯托克斯那个部分。具体推导不谈，反正直到散度的重要性在于，可用表征空间各点矢量场发散的强弱程度

散度的计算如下，$\bold{v} = \{f(x, y, z), g(x, y, z), h(x, y, z)\}$：

\[ \text{div}(\bold{v}) = \frac{\partial f(x, y, z)}{\partial x} + \frac{\partial g(x, y, z)}{\partial y} + \frac{\partial h(x, y, z)}{\partial z} \]

回看 PM 扩散方程的公式，我们先忽略 $c(u)$ 的定义，假设它都是 1；而另一个项目 $\nabla u$ 是梯度，即 $\nabla u = (\partial u / \partial x, \partial u / \partial y)$。

很多图像处理的文章讲到这里就直接给如下公式，然后说这就是拉普拉斯，非常让人迷惑：

\[ \text{div}(\nabla u) = \frac{\partial^{2} u}{\partial x^2} + \frac{\partial^{2} u}{\partial y^2} \]

下面详细讲一下为什么这样，为了方便，把 $\partial u / \partial x$ 记作 $H$，$\partial u / \partial y$ 记作 $V$。

首先明确一点：$H$ 和 $V$ 也都是函数，就像原图 $u(x, y)$ 一样，即他们也是 $H(x, y)$ 和 $V(x, y)$。带入到散度计算公式，得到：
$$
\text{div}(\nabla u)
= \text{div}((H, V))
= \frac{\partial H(x, y)}{\partial x} + \frac{\partial V(x, y)}{\partial y}
$$

现在来看图像中某个像素点 $(p, q)$，它的 $\text{div}(\nabla u)$ 是多少。首先假设做梯度是当前减去左边，当前减去上边：
$$
H(p, q) = u(p, q) - u(p, q-1)
$$
$$
V(p, q) = u(p, q) - u(p-1, q)
$$

那么 $H$ 做 $x$ 方向梯度，得到：
$$
\begin{aligned}
\frac{\partial H}{\partial x}(p, q)
&= H(p, q) - H(p, q-1) \
&= \bigg(u(p, q) - u(p, q-1)\bigg) - \bigg(u(p, q-1) - u(p, q-2)\bigg) \
&= u(p, q) + u(p, q-2) - 2u(p, q-1) \
\end{aligned}
$$

发现结果偏移了一个像素，所以为了方便，规定对 $H$ 做 $x$ 方向梯度时，$H(p, q)$ 是当前减去右边，当前减去下边：
$$
\begin{aligned}
\frac{\partial H}{\partial x}(p, q)
&= H(p, q) - H(p, q+1) \
&= \bigg(u(p, q) - u(p, q-1)\bigg) - \bigg(u(p, q+1) - u(p, q)\bigg) \
&= 2u(p, q) - u(p, q-1) - u(p, q+1) \
\end{aligned}
$$

可以看到正好是中间乘二再减去左边和右边，所以这就是为什么说图像的散度就是拉普拉斯。

因此，可以写为：
$$
\text{div}(\nabla u) = \nabla_{E}(u) + \nabla_{N}(u) + \nabla_{W}(u) + \nabla_{S}(u)
$$
其中 $\nabla_{E}(u) = u(x,y+1)-u(x,y)$，$\nabla_{W}(u) = u(x,y-1)-u(x,y)$，另外两个同理。

扩散传导函数

这里的 $c(x)$ 有不同的类型，常见的有：

\[ c1(x) = \text{exp}\left(-(x / k)^2\right) \]

\[ c2(x) = \frac{1}{1+(x / k)^2} \]

\[ c3(x) = \left\{ \begin{array}{ll} 1 & \text{ if } (x)^2 = 0 \\ 1 - \text{exp}\left(-\frac{3.315}{(x/k)^8}\right) & \text{ if } (x)^2 > 0 \end{array} \right\} \]

这里的 $k$ 表示控制扩散的对比因子，值越大越平滑，越不易保留边缘。具体的函数一般常用前两个。

综合

综上，公式如下，从公式可以发现 $c(u)$ 也分了方向，即不同方向的系数不同：

\[ u(t+1) = u(t) + \lambda * \sum^{E,W,N,S}_{d}\bigg(c(\nabla_d(u)) * \nabla_d(u)\bigg) \]

其中 $c(x)$ 参考上面的扩散传导函数，$\nabla_d(u)$ 是差分算子，$d$ 是方向。

各向异性滤波记录

基本思想

公式解释

迭代

散度计算

扩散传导函数

综合

Comments