用 Python 通过马尔可夫随机场（MRF）与 Ising Model 进行二值图降噪

前言

这个降噪的模型来自 Christopher M. Bishop 的 Pattern Recognition And Machine Learning （就是神书 PRML……），问题是如何对一个添加了一定椒盐噪声（Salt-and-pepper Noise）（假设噪声比例不超过 10%）的二值图（Binary Image）去噪。

原图 -> 添加 10% 椒盐噪声的图：

放在 github 上的可运行完整代码：https://github.com/joyeecheung/simulated-annealing-denoising

建模

下文中的数学表示：

• y_i：噪声图中的像素
• x_i：原图中的像素，对应噪声图中的 y_i

既然噪声图是从原图添加噪声而来，我们拥有了先验知识 1：

y_i和x_i有很强的联系。

一般图片里，每个像素和与它相邻的像素值应当较为接近，比如上图中的黑色笔画和白色负空间，除了边缘以外，黑色的像素周围都是黑色像素，白色像素的周围都是白色像素（连成一片））。这样我们就得到了先验知识 2：

x_i和与它相邻的其他像素也存在较强的联系

如果我们狠一点，假设原图像素只与它的直接相邻像素有联系（即具备条件独立性质），我们就可以得到一个具备局部马尔可夫性质（Local Markov property）的图模型：

在这样一个图模型里，我们有两种团（clique）:

1. {x_i, y_i}，即原图像素与噪声图像素对
2. {x_i, x_j}，其中 x_j 表示与 x_i 相邻的像素

这两种团合并起来，得到的 {x_i, y_i, x_j} 显然是一个最大团（Maximal Clique），此时我们可以利用它来对这个马尔可夫随机场进行 factorization，即求得其联合概率分布关于最大团 x_C = {x_i, y_i, x_j} 的函数。

其中 Z 为 partition function，是 p(x) 的归一化常数（normalization constant），求法参见 PRML 8.3.2。因为与我们的实现不相关，这里不赘述。

ψ_C(x_C) 即所谓的 potential function，为了方便我们通常只求它的对数形式 E(x_C)（按照其物理意义称为 energy function）

关于 factorization 的过程和推导可以参见 PRML 8.3.2，这里我们需要做的是定义一个 energy function。在降噪的过程中 energy 越低，联合概率 P(X=x) （降噪过的图像与原图一致的概率）就越大。

因为我们需要处理的是二值图，首先我们定义 x_i ∈ {-1, +1}，假设这里白色为1，黑色为-1。

对于原图像素与噪声图像素构成的团 {x_i, y_i}，我们定义一项 −ηx_iy_i，其中 η 为一个非负的权重。当两者同号时，这项的值为−η，异号时为 η。这样，当噪声图像素与原图像素较为接近时，这一项能够降低 energy（因为为负）。

对于噪声图中相邻像素构成的团 {x_i, x_j}，我们定义一项 −βx_ix_j，其中 β 为一个非负的权重。这样，当处理过的图像里相邻的像素较为接近时，这一项能够降低 energy（因为为负）。

最后，我们再定义一项 hx_i，使处理过的图像整体偏向某一个像素值。

对图像中的每一个像素，将这三项加起来，就得到我们的 energy function：

对应联合概率

显然 energy 越低，降噪过的图像与原图一致的概率越高。（注意因为我们这里求的 E 已经对整个矩阵求和，即对应 potential function 的积，所以计算联合概率分布的时候不需要再求积）

使用 Python 实现这个 energy function 时，我们可以使用一个 closure 来实现一个 function factory，通过传递beta（β），eta（η）和 h 参数，生成对应的 energy function。此外为了方便，我们假设传入的x和y不是一维向量，而是对应图像的二维矩阵（注意是np.ndarray而不是nd.matrix，前者的*才是array multiplication即逐个元素相乘，后者的*是矩阵乘法）。

import numpy as np

def E_generator(beta, eta, h):
    """Generate energy function E.

    Usage: E = E_generator(beta, eta, h)
    Formula:
        E = h * \sum{x_i} - beta * \sum{x_i x_j} - eta * \sum{x_i y_i}
    """
    def E(x, y):
        """Calculate energy for matrices x, y.

        Note: the computation is not localized, so this is quite expensive.