Bayes，HMM，MRF & Gibbs Distribution在图像降噪中的应用

1 问题描述

假设原始图像${\mathbf{X}}_{m\times n}$，第$i$行，第$j$列的灰度为$x_{i,j}$，经过高斯噪声后
$$\mathbf{Y}=\mathbf{X}+\mathbf{N}$$ 其中$\mathbf{N}$是$m\times n$噪声矩阵，第$i$行，第$j$列的灰度为$N_{i,j}$，$N_{ij}\sim \mathcal{N}(0,{\sigma }^2)$。例如，原始图像$\mathbf{I}$为

添加噪声后得到观测$\mathbf{Y}$

问题是，如何利用像素之间的相互关系，基于$\mathbf{Y}$ 恢复 $\mathbf{X}$

2 数学模型

2.1Bayes推断模型

将降噪问题建模为基于观测$\mathbf{Y}$推断$\mathbf{X}$问题，基于最大后验概率准则（MAP：Maximum A Posteriori)：
$$\hat{\mathbf{X}}=\arg \underset{\mathbf{X}}{\max }\ell =p(\mathbf{X}|\mathbf{Y})$$根据Bayes公式：
$$p(\mathbf{X}|\mathbf{Y})=\frac{p(\mathbf{Y}|\mathbf{X})p(\mathbf{X})}{p(\mathbf{Y})}$$
优化问题等价于：
$$\hat{\mathbf{X}}=\arg \underset{\mathbf{X}}{\max }{\ell }^{\prime }=p(\mathbf{Y}|\mathbf{X})p(\mathbf{X})$$其中似然函数$p(\mathbf{Y}|\mathbf{X})$由噪声模型给出(假设各像素的噪声相互独立)：
$$p(\mathbf{Y}|\mathbf{X})=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\sum _{i=1}^m\sum _{j=1}^n\frac{(y_{i,j}-x_{i,j}{)}^2}{2{\sigma }^2}\right\}$$问题的主要困难有两处：

1. 如何对$\mathbf{X}$的先验分布$p(\mathbf{X})$进行准确建模？
2. $\mathbf{X}$的搜索空间非常庞大，如何高效求解上述优化模型？

2.2 HMM模型

先验分布$p(\mathbf{X})$表示在观测前，对于图像$\mathbf{X}$的了解程度。根据概率链式法则计算先验概率$\mathbf{X}$：
$$p(\mathbf{X})=p(x_n|x_{n-1},x_{n-2},...,x_1)\cdots p(x_3|x_2,x_1)p(x_2|x_1)p(x_1)$$上述展开是很困难的，为了简化模型，可以将概率模型进行简化，认为$\mathbf{X}$的各个分量按顺序形成一条马尔科夫链：

此时有
$$p(\mathbf{X})=p(x_n|x_{n-1})\cdots p(x_3|x_2)p(x_2|x_1)p(x_1)$$ 这就是隐马尔科夫模型 HMM(HMM: Hiden Markov Model)所采用的隐状态概率模型。

HMM 是Bayes推断的一种特例。HMM模型对Bayes推断中的$p(\mathbf{X})$进行简化假设。HMM假设隐状态$\mathbf{X}$的各个随机变量$x_i$是条马尔科夫链，即当前状态$x_i$只与上一个状态$x_{i-1}$有关，而与如何到达上一个状态的过程无关。
HMM的目标是基于观测$\mathbf{Y}$推断$\mathbf{X}$，先验知识是观测模型$\mathbf{Y}=\mathbf{X}+\mathbf{N}$，以及状态转移概率模型$p(x_i|x_{i-1})$，通常HMM采用Bayes方法，即MAP构建目标函数，采用Viterbbi算法求解。Viterbbi算法的基本原理是动态规划（DP: Dynamic Programming）方法。
在语音识别中，常用HMM模型，将音素（phoneme）作为$x_i$，构建相邻两个音素之间的概率关系，构成一条马尔科夫链，运用HMM的方法估计$\mathbf{X}$。
为了进一步贴近事实，也有部分研究对HMM进行了扩展，构建高阶HMM，例如三音素HMM模型，即认为当前音素与前三个因素的状态有关，构建三阶HMM，后面做法与单阶HMM相同。
更多HMM内容参考这里。

2.3 CRF模型

HMM只能表示链状概率转移关系，而有些场景下，链状关系无法准确描述不同随机变量的概率依存关系，这就需要构建MRF（Markov Random Feild）模型了。MRF是链状Markov过程的扩展。MRF中常用农场举例说明。

农场的一个点上种植的作物这个随机变量$x_i$，仅与临近的地种植的作物品种的随机变量$x_j$有关，而且大概率与临近种植作物相同，例如
$$p(x_i=莴苣|x_j=莴苣)=0.8,p(x_i=非莴苣|p_j=莴苣)=0.2)$$
如果每个隐变量还有一个观测$\mathbf{Y}$，我们需要根据观测推断隐变量$\mathbf{X}$的概率分布，即估计后验概率$p(\mathbf{X}|\mathbf{Y})$，则构成条件随机场（CRF：Conditional Random Field）
用概率图表示就是：

图中一个蓝色的圈$x_{ij}$，指的是一个像素的真实状态，对应绿色的圈是该真实状态的观测灰度值$y_{i,j}$。每个蓝色圈对应的像素的真实灰度值与邻近的4个像素的颜色相关，并且被观测到绿色的结果$y_{ij}$

区别于Bayes网络，MRF是一个无向图，适合于解释没有明确因果关系的网络关系。
针对图像而言，先验知识是图像大概率是连续的，光滑的，即图像像素的灰度大概率与临近像素的灰度接近。为了简单，可以假设每个像素的灰度只与4邻域或者8邻域像素的灰度有关，不妨设：
$$p(x_{\bar{i},\bar{j}}|x_{i,j})=\frac{1}{\sqrt{2\pi }{\sigma }_x}\exp \left\{\frac{-(x_{\bar{i},\bar{j}}-x_{i,j}{)}^2}{2{\sigma }_x^2}\right\}$$其中$x_{\bar{i},\bar{j}}$表示第$i$行第$j$列的某个相邻像素的灰度值，例如4邻域中， ( i ˉ , j ˉ ) ∈ { ( i − 1 , j ) , ( i + 1 , j ) , ( i , j − 1 ) , ( i , j + 1 ) } (\bar{i},\bar{j})\in\{(i-1,j),(i+1,j),(i,j-1),(i,j+1)\} (iˉ,jˉ​)∈{(i−1,j),(i+1,j),(i,j−1),(i,j+1)}。${\sigma }_v$是相邻像素灰度差异的标准差，该数值可以通过统计大量自然图像估计得到。
基于上述模型，可以得到一个局部的MRF：

该局部模型的Bayes推断模型为：
$$\begin{gathered} p(y_{i,j}|x_{i,j})=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{(y_{i,j}-x_{i,j}{)}^2}{2{\sigma }^2}\right\} \\ p(x_{\bar{i},\bar{j}}|x_{i,j})=\frac{1}{\sqrt{2\pi }{\sigma }_v}\exp \left\{-\frac{(x_{\bar{i},\bar{j}}-x_{i,j}{)}^2}{2{\sigma }_v^2}\right\} \end{gathered}$$其中前式为观测模型，后式为状态关系模型。
若用Bayes推断的方式估计最优$\mathbf{X}$，需要计算$p(\mathbf{Y}|\mathbf{X})$和$p(\mathbf{X})$。其中，假设观测噪声独立，则有：
$$p(\mathbf{Y}|\mathbf{X})=\sum _{i=1}^m\sum _{j=1}^{n}p(y_{i,j}|x_{i,j})$$现在新的问题是需要根据局部条件概率$p(x_{\bar{i},\bar{j}}|x_{i,j})$，估计联合概率分布$p(\mathbf{X})$。也就是说，需要基于概率图中每条边的概率转移关系，推断联合概率分布。传统方法是基于链式法则展开，但展开式会非常复杂。这里就需要引入Gibbs分布的概念。

2.4 Gibbs分布

MRF研究的两个重要问题是：

1. 如何根据概率图中边的概率关系$p(x_i|x_j)$推断所有随机变量的联合概率$p(\mathbf{X})$；
2. 如何根据概率图的概率关系生成服从$p(\mathbf{X})$分布的随机样本，主要用于马尔科夫蒙特卡洛（MCMC: Markov Chain Monte Carlo)，该问题可以使用metropolis-hasting算法（低维随机数生成有效）和它的特例Gibbs采样算法（主要针对高维随机数生成）解决。

针对第二种方法，在后续文中讨论，这里主要利用MRF研究的第一个问题解决图像降噪问题。
Hammersley Clifford Theorem给出了马尔科夫随机场的联合概率$p(\mathbf{X})$的计算方法，即MRF的联合概率$p(\mathbf{X})$和Gibbs分布是一致的。也就是说，可以使用Gibbs分布估计$p(\mathbf{X})$。

Hammersley Clifford Theorem指出：

1. Gibbs分布一定满足由node separation导致的条件独立性
2. 马尔科夫随机场的概率分布一定可以表示成最大团(Maximum Clique)上的非负函数乘积形式
   证明方法可以看这里

概率图中的”团”指的是概率图中的一个全连接子图，即“团”内的任意两点都有边连接。“最大团”指的是不被任何其他“团”包含的团。在图像降噪问题的例子中，没有两个以上的全连接子图，最大团是包含两个节点一条边的团，所以图像降噪例子的各个像素灰度的联合分布可以表示成每条边上一个非负函数乘积的形式：
$$p(\mathbf{X}|\mathbf{Y})=\frac{1}{z}\prod _{i=1}^m\prod _{j=1}^m\exp F_{i,j}(\mathbf{X},\mathbf{Y}))$$其中$z$是归一化系数，$F_{i,j}(\mathbf{X})$表示与像素$x_{i,j}$有关的势函数，它包含两部分，一部分是像素$x_{i,j}$与观测$y_{i,j}$的势函数$F_o(x_{i,j},y_{i,j})$，另一部分是像素$x_{i,j}$与临近像素$x_{i-1,j}$和$x_{i,j-1}$的势函数$F_g(x_{i,j},x_{i-1,j},x_{i,j-1})$，并且
$$F_{i,j}(\mathbf{X})=F_o(x_{i,j},y_{i,j})F_g(x_{i,j},x_{i-1,j},x_{i,j-1})$$为了避免重复计算，只计算每个像素与左侧和下侧像素的势函数$F_g(x_{i,j},x_{-i,j})$，为了避免第一行和第一列计算异常，定义$x_{i,j}=0$,if $i=0$ or $j=0$）。
在高斯噪声假设和相邻像素灰度差异服从高斯分布的假设下，定义势函数为：
$$\begin{gathered} F_o(x_{i,j},y_{i,j})\triangleq -\frac{1}{2{\sigma }^2}(x_{i,j}-y_{i,j}{)}^2 \\ {F}_g(x_{i,j},x_{i-1,j},x_{i,j-1})\triangleq -\frac{1}{2{\sigma }_v^2}(x_{i,j}-x_{i-1,j}{)}^2+(x_{i,j}-x_{i,j-1}{)}^2 \end{gathered}$$
根据最大后验概率（MAP）准则，对后验概率函数取对数，去除常数项后：
$$\hat{\mathbf{X}}=\arg \min {\ell }^{\prime \prime }=\sum _{i=1}^m\sum _{j=1}^n\left\{{\lambda }_o(x_{i,j}-y_{i,j}{)}^2+{\lambda }_g\left[(x_{i,j}-x_{i-1,j}{)}^2+(x_{i,j}-x_{i,j-1}{)}^2\right]\right\}$$其中${\lambda }_o=\frac{1}{2{\sigma }^2}$表示观测的权重，${\lambda }_g=\frac{1}{2{\sigma }_v^2}$表示图像先验知识的权重，$x_{0,j}\triangleq x_{1,j},x_{i,0}\triangleq x_{i,1}$。事实上，对比Bayes推断框架，前者对应观测似然函数$p(\mathbf{Y}|\mathbf{X})$，后者对应先验知识似然函数$p(\mathbf{X})$。

经过上述过程，构建起了推断$\mathbf{X}$的数学模型。在求解该优化问题时，就需要遍历所有的$\mathbf{X}$，找出最大后验的解，此时搜索规模是$256^{m\times n}$，显然是无法接受的，需要设计有效的搜索算法。

3 快速算法

优化目标是搜索最佳的$x_{i,j}$，使两类距离（观测与估计距离&相邻像素之间距离）最小。若能每步迭代都减小一点距离，并且有个好的初值，就能以比较大的概率收敛。考虑以下迭代方法：
$${\hat{x}}_{i,j}^{(t+1)}=\arg \underset{x_{i,j}^{(t+1)}}{\min }{\ell }_{i,j}^{(t+1)}={\lambda }_o(x_{i,j}-y_{i,j}{)}^2+{\lambda }_g\left[(x_{i,j}-x_{i-1,j}^{(t)}{)}^2+(x_{i,j}-x_{i,j-1}^{(t)}{)}^2\right]\text{\hspace{1em}(1)}$$ 直观意义就是优化第$i$行第$j$列的灰度值$x_{i,j}$，使其与观测$y_{i,j}$和临近像素灰度值$x_{i-1,j},x_{i,j-1}$尽可能接近。基于Gibbs抽样的相关结论，不难证明，该迭代过程是收敛于MAP估计的（具体证明方法在后续文章给出）。
一个好的初值能够加速收敛过程，注意到真实图像$\mathbf{X}$应该接近观测$\mathbf{Y}$，所以使用$\mathbf{Y}$作为$\mathbf{X}$的初值是合理的。
综上得到算法流程：

其中逐像素更新模块为求解上述优化模型(1)，迭代终止条件为图像灰度值不再发生变化。

4 代码与结果分析

代码如下，代码部分参考了这里：

% close all;
clear all;
figure;
I=rgb2gray(imread('lena.jpg'));
imwrite(I,'lena_gray.jpg');
I1=double(I);
subplot(2,2,1)
imagesc(I1);
title('原图像');
colormap('gray')
J1 = double(imnoise(I,'gaussian',0,0.001));
subplot(2,2,2)
imagesc(J1);
title('噪声图')
colormap('gray')
imwrite(uint8(J1),'lena_noised.jpg');

Y = J1;
[m,n]=size(Y);
X = Y;
h=0;beta=0.05;eta=1;
z = [-1,0,1];
while 1
    tot = 0;
    for i=2:1:m-1
        for j=2:1:n-1
            temp = X(i,j);
            sq = X(i-1:i+1,j-1:j+1);
            sq1 = sq;
            sq1(2,2) = sq(2,2)-1;   %根据定义计算势函数
            E1 = beta*((sq1(2,2)-sq1(1,2))^2+(sq1(2,2)-sq1(2,1))^2+(sq1(2,2)-sq1(2,3))^2+(sq1(2,2)-sq1(3,2))^2)+eta*(sq1(2,2)-Y(i,j))^2/9;
            sq2 = sq;
            sq2(2,2) = sq(2,2);   %根据定义计算势函数
            E2 = beta*((sq2(2,2)-sq2(1,2))^2+(sq2(2,2)-sq2(2,1))^2+(sq2(2,2)-sq1(2,3))^2+(sq2(2,2)-sq2(3,2))^2)+eta*(sq2(2,2)-Y(i,j))^2/9;
            sq3 = sq;
            sq3(2,2) = sq(2,2)+1;   %根据定义计算势函数
            E3 = beta*((sq3(2,2)-sq3(1,2))^2+(sq3(2,2)-sq3(2,1))^2+(sq3(2,2)-sq3(2,3))^2+(sq3(2,2)-sq3(3,2))^2)+eta*(sq3(2,2)-Y(i,j))^2/9;
            [~,index] = min([E1,E2,E3]);
            X(i,j) = X(i,j) + z(index);
            if temp~=X(i,j)
            tot=tot+1;
            end
        end
    end
    if tot<1
        break;
    end
end
 
J2=X;
subplot(2,2,3)
imagesc(J2);
colormap('gray')
title('mrf降噪结果')

subplot(2,2,4)
imagesc(J2-J1);
colormap('gray')
title('差异')
imwrite(uint8(J2),'lena_denoised.jpg');

结果

原始图

加噪声图

降噪图

从结果看出，噪声得到一定程度的遏制，但是图像锐度发生下降，这是因为代价函数总是趋向于让图像变得光滑导致的。如果把更准确的先验知识应用代价函数中，则可以获得更理想的效果。
例如为了避免图像变化剧烈的正常边缘被模糊，将先验知识中相邻像素的势函数修改为
$$F_g(x_{i,j},x_{i-1,j},x_{i,j-1})\triangleq -\frac{1}{2{\sigma }_v^k}|(x_{i,j}-x_{i-1,j}){|}^k+|(x_{i,j}-x_{i,j-1}){|}^k$$
$k$越小，则表示越忽视过大差异的势，若边缘模糊，则说明代价函数过多关注了相邻像素的光滑性，需要将k减小，以提高图像的锐度，例如当$k=1.5$时，图像锐度得到了提高：

上述三张图分别时含噪图，$k=2$ 和$k=1.5$，可以看出$k$的下降有利于保护原图中高反差的边缘（例如胳膊和背景，帽子亮部和背景），但一些低反差的细节可能会损失（例如发丝细节）。实际中可以通过机器学习的方式获得最佳的参数估计。