基于分塊排序重采樣PCA的泊松降噪算法

【作 者】郭哲，趙文釗，秦斌杰

上海交通大學生物醫(yī)學工程學院，上海市，200240

基于分塊排序重采樣PCA的泊松降噪算法

【作 者】郭哲，趙文釗，秦斌杰

上海交通大學生物醫(yī)學工程學院，上海市，200240

泊松噪聲在低光子計數成像中較為常見，尤其是在微光成像、天文、核醫(yī)學領域，但由于建模處理信號相關噪聲的困難性，對于低光子計數(及小尺寸圖像)往往會存在小樣本問題以及圖像區(qū)塊間特征自相似性不足的問題，使得當今許多優(yōu)秀的降噪算法還不能達到好的降噪效果。該文提出了一種能同時解決這兩種問題的泊松降噪算法。首先，我們對圖像塊進行分塊排序，用重采樣法對鄰近非局部塊進行分塊重采樣；然后選取與原始圖像塊近似度高的前5個采樣向量并應用基于指數分布族的PCA框架對其進行處理；最后依據與原始圖像近似度權重將處理結果合成降噪后的圖像。實驗結果顯示我們的方法對小樣本量低光子計數圖像的泊松降噪有很好的表現。

泊松噪聲；低光子計數問題；小樣本問題；分塊排序；PCA；重采樣

0 引言

在微弱光成像、X射線光譜成像、熒光顯微成像和天文學成像時，在有限的曝光時間內，只有有限數量的光子數被采集。這些低光子計數光學圖像具有極低的信噪比[1]。由于量子力學的不確定性原則，由樣本中輻射出的光子的波動是隨機變化的，所以傳感器上不同像素點檢測到的光子數可以看作是遵循泊松分布的獨立事件，低光子計數圖像往往會受到泊松噪聲污染。由于建模處理信號相關噪聲的困難性，傳統(tǒng)的降噪算法往往忽略考慮泊松噪聲特性，只是對加性高斯白噪聲模型進行降噪處理。

目前，三維塊匹配(BM3D)[2]方法被認為是當今對降噪最為有效的算法，通過將二維圖像塊變成三維數據序列來增加信號表示的稀疏性，結合協(xié)同維納濾波來有效地降低噪聲。盡管這個算法對于信號獨立的高斯噪聲表現良好，但其并沒有考慮現實中的低光子計數圖像更應該遵循泊松噪聲模型[3]。由于泊松噪聲是信號相關的，我們就無法估計所期望的噪聲方差常量，這就為設計降噪算法帶來了困難。為了克服這些困難，方差穩(wěn)定變換(VST)諸如Anscombe 變換[4]就被引入來處理泊松噪聲。Donoho利用Anscombe VST進行了泊松降噪[5]。由于Anscombe 變換是非線性變換，反變換帶來的偏離誤差是不可避免的。為了克服這種偏離，Makitalo 和 Foi提出了廣義Anscombe變換(GAT)[6-7]。另外一類泊松噪聲的降噪算法使用了Haar變換，結合假設檢驗[8]和貝葉斯框架[9-10]得到了很大的應用。Kolaczyk開發(fā)了收縮法修正基于任意小波變換的硬/軟閾值來處理泊松噪聲[11]。

最近又出現了一種針對泊松噪聲的稀疏正則化約束的凸優(yōu)化算法[12]，Salmon等用針對泊松噪聲的非局部主成分分析法改進了這一算法[13]。Giryes和Elad提出的泊松降噪方法，采用了Salmon等提出的方法再結合了字典學習的稀疏表示[14]。Elad 和 Anaron[15]提出使用正交匹配(OMP)進行稀疏編碼學習的過完備字典完成圖像的降噪重建。但是這些字典學習方法沒有考慮低光子計數圖像的小樣本問題。小樣本問題對基于字典學習的泊松降噪帶來極大挑戰(zhàn)。字典學習方法總是假定圖像數據特定的顯著特征可以通過大量樣本的訓練而可以被一個合適的字典表達[16]。但小樣本問題使得每個圖像塊的多種顯著特征被各不相同的泊松分布所表示，通過有限計數的小樣本訓練表達這些特征變得十分困難。

因此，本文中，我們提出了一個簡單但有效的基于泊松噪聲模型的降噪算法，解決存在小樣本問題的低光子計數圖像降噪問題。我們認為泊松計數采集過程的異方差性決定了每個像素位置極其有限的光子計數基本上和周圍領域像素分屬于不同的獨立統(tǒng)計分布，因此決定了低光子計數過程存在了小樣本問題。我們對提取出來的各圖像塊進行排序，然后用自助重采樣法對這些圖像塊進行重采樣；然后選取與原始圖像塊近似度高的前5個采樣向量并應用基于指數分布族的PCA框架對其進行處理；最后依據與原始圖像近似度權重將處理結果合成降噪后的圖像。我們的實驗數據表明，我們的算法比當今的泊松降噪算法還有較好的表現。

1 泊松噪聲模型

基于光子計數的光譜圖像數據是典型的泊松分布數據。假設zi（i=1,..., N，N是掃描點陣的總數目）是在圖像采集設備上觀測到的像素強度值，zi服從函數的統(tǒng)計分布，其中λ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像素的圖像塊進行算法降噪性能的比較。我們之所以考慮小尺寸圖像是由于它們更存在了小樣本和區(qū)塊間圖像特征自相似性不足的問題。我們對這些圖像加入泊松噪聲(其中噪聲峰值分別為peak=0.1以及peak=1)，得到噪聲圖像。然后利用本文方法對圖像進行降噪處理。為了比較我們提出算法的性能優(yōu)劣，實驗中還采用了結合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中看出，本方法綜合地優(yōu)于其他兩種算法，特別是在對小尺寸圖像進行降噪處理，更能體現本文算法的優(yōu)勢。這是因為我們的算法在重采樣操作以后，更能充分捕捉圖像本身的統(tǒng)計特性和所表征的圖像結果特征。

表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 總結

我們的算法首先使用塊排序操作，將圖像中相似的塊排布在一起，然后利用重采樣的方法來擴大數據量，找到更多的和原始圖像塊更為相似的圖像塊矩陣；緊接著通過結合隨機距離的思想和主成分分析算法，能夠有效地去除圖像中存在的泊松噪聲，達到較好的降噪性能。實驗結果表明，本文提出的算法較其他算法而言具有更好的降噪性能，尤其是當圖像的尺寸偏小時，我們的算法更具有明顯的優(yōu)勢。我們的算法仍有一些不足：如算法對圖像細節(jié)修復的不夠好，所以我們需要進一步地優(yōu)化算法，以提高算法的降噪性能；目前，我們的實驗數據還不是很充分，后續(xù)還要大量的小尺寸數據進行實驗。

[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）；上海交通大學醫(yī)工交叉基金面上項目(YG2014MS29)

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