【图像处理】基于贝叶斯的Lucy

研究空间域滤波方法是图像处理领域的重要内容。相比频域滤波，空间域滤波直接在图像空间上进行操作，除了传统的滤波外，还能加入各种直观的空间域操作，可扩展性和可理解性都比较强。在图像复原领域，当噪声较为复杂时，频域滤波方法因为难以计算得到噪声的频域模型，不能合理地去除噪声干扰，空间域方法就容易多了，因此空间域滤波方法占有很重要的地位。

图像复原问题的重要步骤是建立合理的图像退化模型。假设退化图像，即观测图像为H，原图为W，点扩散函数（可以理解为模糊卷积核）为S，噪声为N，*为卷积符号，则退化模型为 H = W ∗ S + N H=W*S+N H=W∗S+N，所求为W。若假设图像误差满足高斯概率分布，那么可以用最小二乘思想来优化解决上述退化模型，模型为 a r g m i n 1 2 ( W ∗ S − H ) 2 argmin \frac{1}{2}(W*S-H)^2 argmin21​(W∗S−H)2，现有方法比如维纳滤波、最小二乘空间域滤波都适用。不过如果要把退化模型推广到任意概率分布，就要寻找另外的更鲁棒的方法框架。

Lucy-Richardson图像复原方法（LR方法或者RL方法）是一种基于贝叶斯思想的空间域上的图像复原方法。该算法从贝叶斯理论出发推导了图像迭代复原的基本框架。本文根据Richardson W H. 1972. 给出一个简单的例子。

我们知道光学图像中，由于成像时光线有散射等情况，导致图像中某个像素的像素值，由本身像素点的强度和邻近像素点的光线干扰组成。设原始图像第i个像素点为 W i W_i Wi​，观测图像第k个像素点为 H k H_k Hk​， W i W_i Wi​对 H k H_k Hk​的像素值贡献可以表示为 S i , k S_{i,k} Si,k​。

如果把图像归一化，像素值当作概率值，那么整个图像辨成一个概率图， W i W_i Wi​等价于 P ( W i ) P(W_i) P(Wi​)， S i , k S_{i,k} Si,k​等价于 P ( H k ∣ W i ) P(H_k|W_i) P(Hk​∣Wi​)，那么反过来已知观测图像，潜在的原始图像分布概率为 P ( W i ∣ H k ) = P ( H k ∣ W i ) P ( W i ) ∑ j P ( H k ∣ W j ) P ( W j ) (1) P(W_i|H_k)=\frac{P(H_k|W_i)P(W_i)}{\sum_j{P(H_k|W_j)P(W_j)}} \tag{1} P(Wi​∣Hk​)=∑j​P(Hk​∣Wj​)P(Wj​)P(Hk​∣Wi​)P(Wi​)​(1)，其中j代表点扩散函数S的作用范围； P ( H k ∣ W i ) P(H_k|W_i) P(Hk​∣Wi​) 表示点扩散函数，可写作 P ( S i , k ) P(S_{i,k}) P(Si,k​)。注意这里，归一化图像等效为概率图，并不是说某个像素值出现的概率服从某个概率分布。实际上是用贝叶斯原理等效解决图像空间域的问题。

设原图为 P ( W i ) P(W_i) P(Wi​)，写出一个取巧的联合分布公式 P ( W i ) = ∑ k P ( W i , H k ) = ∑ k P ( W i ∣ H k ) P ( H k ) (2) P(W_i)=\sum_k{P(W_i,H_k)}=\sum_k{P(W_i|H_k)P(H_k)} \tag{2} P(Wi​)=k∑​P(Wi​,Hk​)=k∑​P(Wi​∣Hk​)P(Hk​)(2)右侧是全概率公式转变而来的。联合公式(1)(2)可得： P ( W i ) = ∑ k P ( H k ∣ W i ) P ( W i ) P ( H k ) ∑ j P ( H k ∣ W j ) P ( W j ) (3) P(W_i) = \sum_k{\frac{P(H_k|W_i)P(W_i)P(H_k)}{\sum_j{P(H_k|W_j)P(W_j)}}} \tag{3} P(Wi​)=k∑​∑j​P(Hk​∣Wj​)P(Wj​)P(Hk​∣Wi​)P(Wi​)P(Hk​)​(3)注意内外求和符号中，求和的对象是j和k。(3)中左右两侧都有 P ( W i ) P(W_i) P(Wi​)，联系牛顿下降法的推导，可以得到一个迭代公式： P r + 1 ( W i ) = ∑ k P ( H k ∣ W i ) P ( H k ) P r ( W i ) ∑ j P ( H k ∣ W j ) P r ( W j ) (4) P^{r+1}(W_i) = \sum_k{\frac{P(H_k|W_i)P(H_k)P^r(W_i)}{\sum_j{P(H_k|W_j)P^r(W_j)}}} \tag{4} Pr+1(Wi​)=k∑​∑j​P(Hk​∣Wj​)Pr(Wj​)P(Hk​∣Wi​)P(Hk​)Pr(Wi​)​(4)

现将(4)应用到真实图像复原场景中去。Richardson W H. 1972. 设定 P ( W i ) = W i / W P(W_i)= W_i/W P(Wi​)=Wi​/W， W W W是归一化系数，对其他项也是同样的。(4)中各项概率值可以直接等效为像素值，其中归一化系数在上下分式中被舍弃： W i r + 1 = ∑ k S i , k H k W i r ∑ j S j , k W j r (5) W_i^{r+1} = \sum_k{\frac{S_{i,k}H_kW^r_i}{\sum_j{S_{j,k}W^r_j}}} \tag{5} Wir+1​=k∑​∑j​Sj,k​Wjr​Si,k​Hk​Wir​​(5)

现将(5)扩展到二维图像复原。首先必须注意到内外求和公式中求和项分别是k和j。(5)中外面的求和项中，i是当前像素点，是当前处理区域中心点；而在内部下方的求和中，j是可变的，k才是中心点。因此图像复原公式为： W r + 1 = W r S f l i p ∗ H S ∗ W r W^{r+1}=W^r\frac{S^{flip}*H}{S*W^r} Wr+1=WrS∗WrSflip∗H​

实际应用中，观测图像、点扩散函数是已知的。本文对lena图像进行了11x11的高斯滤波，并就加入高斯噪声、盐噪声进行了测试。 高斯滤波+高斯噪声 高斯滤波+高斯噪声的复原效果，细节还原还好，振铃效应显著，需要改进的LR算法。

高斯滤波+5%盐噪声 高斯滤波+5%盐噪声的复原结果，效果很差，说明原始LR算法对噪声是很敏感的。 对观测图像（高斯滤波+5%盐噪声）进行3x3中值滤波之后的复原结果，效果比直接用LR复原要好得多。

[1] Richardson W H . Bayesian-Based Iterative Method of Image Restoration*[J]. Journal of the Optical Society of America (1917-1983), 1972, 62(1):55-59.