基于分块排序重采样PCA的泊松降噪算法

【作 者】郭哲，赵文钊，秦斌杰

上海交通大学生物医学工程学院，上海市，200240

基于分块排序重采样PCA的泊松降噪算法

【作 者】郭哲，赵文钊，秦斌杰

上海交通大学生物医学工程学院，上海市，200240

泊松噪声在低光子计数成像中较为常见，尤其是在微光成像、天文、核医学领域，但由于建模处理信号相关噪声的困难性，对于低光子计数(及小尺寸图像)往往会存在小样本问题以及图像区块间特征自相似性不足的问题，使得当今许多优秀的降噪算法还不能达到好的降噪效果。该文提出了一种能同时解决这两种问题的泊松降噪算法。首先，我们对图像块进行分块排序，用重采样法对邻近非局部块进行分块重采样；然后选取与原始图像块近似度高的前5个采样向量并应用基于指数分布族的PCA框架对其进行处理；最后依据与原始图像近似度权重将处理结果合成降噪后的图像。实验结果显示我们的方法对小样本量低光子计数图像的泊松降噪有很好的表现。

泊松噪声；低光子计数问题；小样本问题；分块排序；PCA；重采样

0 引言

在微弱光成像、X射线光谱成像、荧光显微成像和天文学成像时，在有限的曝光时间内，只有有限数量的光子数被采集。这些低光子计数光学图像具有极低的信噪比[1]。由于量子力学的不确定性原则，由样本中辐射出的光子的波动是随机变化的，所以传感器上不同像素点检测到的光子数可以看作是遵循泊松分布的独立事件，低光子计数图像往往会受到泊松噪声污染。由于建模处理信号相关噪声的困难性，传统的降噪算法往往忽略考虑泊松噪声特性，只是对加性高斯白噪声模型进行降噪处理。

目前，三维块匹配(BM3D)[2]方法被认为是当今对降噪最为有效的算法，通过将二维图像块变成三维数据序列来增加信号表示的稀疏性，结合协同维纳滤波来有效地降低噪声。尽管这个算法对于信号独立的高斯噪声表现良好，但其并没有考虑现实中的低光子计数图像更应该遵循泊松噪声模型[3]。由于泊松噪声是信号相关的，我们就无法估计所期望的噪声方差常量，这就为设计降噪算法带来了困难。为了克服这些困难，方差稳定变换(VST)诸如Anscombe 变换[4]就被引入来处理泊松噪声。Donoho利用Anscombe VST进行了泊松降噪[5]。由于Anscombe 变换是非线性变换，反变换带来的偏离误差是不可避免的。为了克服这种偏离，Makitalo 和 Foi提出了广义Anscombe变换(GAT)[6-7]。另外一类泊松噪声的降噪算法使用了Haar变换，结合假设检验[8]和贝叶斯框架[9-10]得到了很大的应用。Kolaczyk开发了收缩法修正基于任意小波变换的硬/软阈值来处理泊松噪声[11]。

最近又出现了一种针对泊松噪声的稀疏正则化约束的凸优化算法[12]，Salmon等用针对泊松噪声的非局部主成分分析法改进了这一算法[13]。Giryes和Elad提出的泊松降噪方法，采用了Salmon等提出的方法再结合了字典学习的稀疏表示[14]。Elad 和 Anaron[15]提出使用正交匹配(OMP)进行稀疏编码学习的过完备字典完成图像的降噪重建。但是这些字典学习方法没有考虑低光子计数图像的小样本问题。小样本问题对基于字典学习的泊松降噪带来极大挑战。字典学习方法总是假定图像数据特定的显著特征可以通过大量样本的训练而可以被一个合适的字典表达[16]。但小样本问题使得每个图像块的多种显著特征被各不相同的泊松分布所表示，通过有限计数的小样本训练表达这些特征变得十分困难。

因此，本文中，我们提出了一个简单但有效的基于泊松噪声模型的降噪算法，解决存在小样本问题的低光子计数图像降噪问题。我们认为泊松计数采集过程的异方差性决定了每个像素位置极其有限的光子计数基本上和周围领域像素分属于不同的独立统计分布，因此决定了低光子计数过程存在了小样本问题。我们对提取出来的各图像块进行排序，然后用自助重采样法对这些图像块进行重采样；然后选取与原始图像块近似度高的前5个采样向量并应用基于指数分布族的PCA框架对其进行处理；最后依据与原始图像近似度权重将处理结果合成降噪后的图像。我们的实验数据表明，我们的算法比当今的泊松降噪算法还有较好的表现。

1 泊松噪声模型

基于光子计数的光谱图像数据是典型的泊松分布数据。假设zi（i=1,..., N，N是扫描点阵的总数目）是在图像采集设备上观测到的像素强度值，zi服从函数的统计分布，其中λi是没有噪声的真实光子计数。泊松变量的均值和方差为

因此，被泊松噪声污染的图像降噪问题就可以等效为从噪声图像中估计无噪图像真实光子计数 λi。

2 基于分块排序重采样PCA的降噪算法

2.1 分块排序算法

对于低光子计数的光谱图像而言，我们先利用图像的分块算法对光谱图像进行分块处理：

其中，x表示的是需要处理的图像，Ri表示的是在某个像素点i处从原始图像x中提取图像块（块的尺寸为×）的矩阵提取操作，xi表示的是将提取出来的块重新排列成一个列向量（尺寸为n），假设图像的尺寸为M×N，则总共可以提取出个不同的图像块。

对图像进行分块以列向量的形式堆叠成矩阵X，然后对这个图像块矩阵进行排序，即寻找一个变换矩阵P使得Xp=PX。通常假定对于干净图像Xp是平滑的，Ram I等[17]基于这一假设提出了分块排序的降噪算法，用总变分作为“平滑性”的测量标准：

最小化‖Xp‖TV需要找到一个最短的路径使得每个点只被访问一次，Ram I采用一种近似的解决方案，具体为：从任意点开始，每一点xj0到下一点xj1的概率为：

到第二近的点xj2的可能性为：

这里ε是预先设定的参数，α被选取适当的值使得p1+p2=1，xj1和xj2取自未被访问的点集。每一点的搜索范围限定在以该点为中心的一定大小的方块中，如果该方块中没有访问的点则在整个图像的未访问点中寻找最近邻点。本文中我们引入这种排序方法进行分块排序。排序完成后用每一点两侧的小块拼接成2N×1的向量yi代替原来的N×1的小块xi，供下一步重采样。

2.2 重采样方法

对提取出的图像块向量进行重采样处理，以获得更多的样本数据。重采样的方法可以用以下公式表示：

2.3 权重计算

我们可以利用NLM算法中权重的思想来对光谱图像进行降噪处理。NLM[17]算法中像素间的权重是通过式(7)来进行计算的：

其中，h是滤波器参数，其控制的是滤波器的延滞操作。当h→0，降噪后的图像更加趋近于输入的噪声图像；当h→∞，则输出图像变得更为模糊。Wi是一个正则化因子，是所有权重的总和Wi=∑j∈Iωij。zηi和zηj分别代表的是把块ηi和ηj展开为列向量的强度值，‖zηi-zηj表示的是利用像素间的欧式距离来进行权重计算。

在本算法中，我们利用随机距离[18]的思想来选取近似图像块。假设两个像素点X和Y分别是服从参数为λX和λY的泊松随机变量，则这两个像素点X和Y之间的距离可以表示为

其中，disXY表示的是像素点X和Y之间的距离。

因此，参照式（8），对于上述通过重采样得到的图像块样本矩阵，我们利用随机距离来计算重采样样本矩阵每一列向量数据与原始被采集的子块向量数据之间的权重关系，即

其中，dis(a, b)用式（8）来表示，Wi表示的是所有权重的总和Wi=∑j∈Iωij，λ(ηi)表示的是原始被采集的向量的期望值，λ(Yjxi)表示的是重采样得到的样本矩阵Yjxi的每一个列向量的期望值，γ为指数的控制因子，设为常数。然后对这些求到的权重数据从大到小地进行排序，获取前M(这里为5)个最大权重所对应的样本矩阵中的列向量，这些列向量就可以被选取为由块xi所提取的近似图像块Dxi，即：

其中，Dxi为选取的近似图像块矩阵，为前M个较大的权重，Yjxi为块xi通过重采样得到的样本矩阵。因此，图像中每个块的近似块经过这样的选取，可以构成一个最终的近似图像块矩阵D，即

2.4 基于指数分布族的PCA算法

我们假定重采样得到的相似图像块仍然是服从泊松分布的。对于指数分布的泊松噪声图像数据而言，基于NLPCA算法[13]的思想，该噪声数据通常可以使用以下公式来进行分解：

因此，我们的目的是要最小化式（14）来得到最终的系数矩阵U*和数据主要成分矩阵V*：

因此，最后得到的降噪图像块矩阵为：

然后我们再利用式(3)的逆操作对求得的图像块矩阵进行还原，得到最终的降噪图像Y。为了达到较好的降噪效果，我们对第一次降噪后的图像结果再次进行相同的处理，经过多次迭代，最终达到理想的效果。

3 实验结果

由于我们的算法是针对小尺寸图像进行处理的，因此，我们256×256的标准Saturn灰度图像中截取出一块50×50像素的图像块进行算法降噪性能的比较。我们之所以考虑小尺寸图像是由于它们更存在了小样本和区块间图像特征自相似性不足的问题。我们对这些图像加入泊松噪声(其中噪声峰值分别为peak=0.1以及peak=1)，得到噪声图像。然后利用本文方法对图像进行降噪处理。为了比较我们提出算法的性能优劣，实验中还采用了结合Anscombe变换[6]的BM3D算法[2]以及NLPCA算法[13]进行降噪处理。图1表示的是不同降噪方法进行降噪所得的结果。图中(a)到(e)分别表示的是无噪图像，噪声图像，NLPCA方法的降噪结果，BM3D算法降噪结果以及我们算法的降噪结果。从直观的视觉评判可以看出，我们的算法较其他两种降噪算法具有更好的降噪性能。

图1 从标准Saturn图像中截取出50×50的图像的降噪处理结果Fig.1 The denoised results for 50×50 seized by standard Saturn image

为了能够客观地衡量算法的降噪性能，表1给出了用峰值信噪比（Peak Signal-to-noise Ratios，PSNR）衡量的性能指标。从表1中看出，本方法综合地优于其他两种算法，特别是在对小尺寸图像进行降噪处理，更能体现本文算法的优势。这是因为我们的算法在重采样操作以后，更能充分捕捉图像本身的统计特性和所表征的图像结果特征。

表1 不同算法降噪结果的PSNR比较（dB）Table 1 The PSNR (dB) results of different algorithms

此外，我们使用真实的老鼠图像来对我们的算法性能进行评估。如图2所示。从图2中我们可以看出在针对小尺寸图像的情况下，相较于其他两种算法，我们可以得到一个比较不错的降噪性能。

图2 截取出的老鼠图像降噪结果Fig.2 The denoised results of seized image of mouse

4 总结

我们的算法首先使用块排序操作，将图像中相似的块排布在一起，然后利用重采样的方法来扩大数据量，找到更多的和原始图像块更为相似的图像块矩阵；紧接着通过结合随机距离的思想和主成分分析算法，能够有效地去除图像中存在的泊松噪声，达到较好的降噪性能。实验结果表明，本文提出的算法较其他算法而言具有更好的降噪性能，尤其是当图像的尺寸偏小时，我们的算法更具有明显的优势。我们的算法仍有一些不足：如算法对图像细节修复的不够好，所以我们需要进一步地优化算法，以提高算法的降噪性能；目前，我们的实验数据还不是很充分，后续还要大量的小尺寸数据进行实验。

[1] Luisier F, Blu T, Unser M. Image denoising in mixed Poisson–Gaussian noise[J]. IEEE Trans Imag Proc, 2011, 20(3): 696-708.

[2] Dabov K , Foi A , Katkovnik V , et al. Image denoising by sparse 3-D transform-domain collaborative fltering[J]. IEEE Trans Imag Proc, 2007, 16(8):2080-2095.

[3] Foi A, Trimeche M, Katkovnik V, et al. Practical Poissonian-Gaussian noise modeling and ftting for single-image raw-data[J]. IEEE Trans Imag Proc,2008, 17(10):1737-1754.

[4] Anscombe FJ, Anscombe FJ. The transformation of Poisson, binomial and negative binomial data[J]. Biometrika, 1957, 35(3-4):246-254.

[5] Donoho DL. Nonlinear wavelet methods for recovery of signals, densities, and spectra from indirect and noisy data[J]. Proc Symp Appl Math, 1970:173-205.

[6] Ma kitalo M, Foi A. Optimal inversion of the Anscombe transformation in low-count Poisson image denoising[J]. IEEE Trans Imag Proc, 2009, 20(1): 26-32.

[7] Ma kitalo M, Foi A. A closed-form approximation of the exact unbiased inverse of the Anscombe variance-stabilizing transformation[J]. IEEE Trans Imag Proc, 2011, 20(9): 2697-2698.

[8] Kolaczyk ED. Nonparametric estimation of intensity maps using Haar wavelets and Poisson noise characteristics[J], Astrophys J, 2000, 534: 490-505.

[9] Lefkimmiatis S, Maragos P, Papandreou G. Bayesian inference on multiscale models for Poisson intensity estimation: Applications to photon-limited image denoising[J]. IEEE Trans Imag Proc, 2009, 18:1724-1741.

[10] Timmermann KE, Nowak RD. Multiscale modeling and estimation of Poisson processes with application to photon-limited imaging[J]. IEEE Trans Inf Theory, 1999, 45(2): 846-862.

[11] Kolaczyk ED. Wavelet shrinkage estimation of certain Poisson intensity signals using corrected thresholds[J]. Sta Sin, 1999, 9: 119-135.

[12] Harmany Z, Marcia R, Willett R. This is SPIRAL-TAP: Sparse poisson intensity reconstruction algorithms–theory and practice[J]. IEEE Trans Imag Proc, 2012, 21(3): 1084-1096.

[13] Salmon J, Harmany Z, Deledalle CA, et al. Poisson noise reduction with non-local PCA[J]. J Math Imag Vis, 2014, 48: 279-294.

[14] Giryes R, Elad M. Sparsity based Poisson denoising with dictionary learning[J]. IEEE Trans Imag Proc, 2014, 23(12): 5057-5069.

[15] Aharon M, Elad M, Bruckstein A K. K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation, IEEE Trans. Signal Processing, 54, 4311-4322[J]. IEEE Trans Sign Proc, 2006, 54(11):4311-4322.

[16] Meng D, Zhao Q, Leung Y, et al. Learning dictionary from signals under global sparsity constraint[J]. Neurocomputing, 2013, 119(16):308-318.

[17] Ram I, Elad M, Cohen I. Image processing using smooth ordering of its patches[J]. IEEE Trans Imag Proc, 2013, 22:2764-2774.

[18] Deledalle CA, Tupin F, Denis L. Poisson NL means: Unsupervised non local means for Poisson noise[C]. IEEE Int Conf Imag Proc, 2010: 801-804.

[19] Bindilatti AA, Mascarenhas NDA. A Nonlocal Poisson Denoising Algorithm Based on Stochastic Distances[J]. IEEE Sign Proc Lett, 2013, 20(11):1010-1013.

Poisson Noise Removal Using Patch-order Resampling PCA Algorithm

【 Writers 】GUO Zhe, ZHAO Wenzhao, QIN Binjie
School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai, 200240

The problem of Poisson denoising is common in various photon-limited imaging applications, especially in low-light imaging, astronomy and nuclear medical applications. Due to the small sample problem and the related insufficient self-similarity between patches of whole image, many denoising algorithms cannot obtain the favorable denoising performance. We propose patch-order resampling PCA algorithm for Poisson noise reduction. Firstly, we use the patchordered operations to sort the extracted image patches and exploit the bootstrap resampling method to resample the different blocks of spectral images to obtain more data matrix of image samples. Then, we select the patches with largest weights corresponding to the vectors of image samples data matrix as the most similar patches. Finally, we use principal component analysis algorithm for processing the image to obtain the fnal denoised image. Experiments results show that the proposed method achieves excellent Poisson noise removal performance in the photon-limited images with small sample problems.

Poisson noise, low photon counts, small sample problem, patch order, PCA, bootstrap resampling

R318.6

A

10.3969/j.issn.1671-7104.2016.06.003

1671-7104(2016)06-0403-04

2016-05-26

国家自然科学基金委面上项目（61271320，60872102）；上海交通大学医工交叉基金面上项目(YG2014MS29)

秦斌杰，E-mail: bjqin@sjtu.edu.cn