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

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

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

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

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

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

泊松噪聲；低光子計(jì)數(shù)問(wèn)題；小樣本問(wèn)題；分塊排序；PCA；重采樣

0 引言

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

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

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

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

1 泊松噪聲模型

基于光子計(jì)數(shù)的光譜圖像數(shù)據(jù)是典型的泊松分布數(shù)據(jù)。假設(shè)zi（i=1,..., N，N是掃描點(diǎn)陣的總數(shù)目）是在圖像采集設(shè)備上觀測(cè)到的像素強(qiáng)度值，zi服從函數(shù)的統(tǒng)計(jì)分布，其中λi是沒(méi)有噪聲的真實(shí)光子計(jì)數(shù)。泊松變量的均值和方差為

因此，被泊松噪聲污染的圖像降噪問(wèn)題就可以等效為從噪聲圖像中估計(jì)無(wú)噪圖像真實(shí)光子計(jì)數(shù) λi。

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

2.1 分塊排序算法

對(duì)于低光子計(jì)數(shù)的光譜圖像而言，我們先利用圖像的分塊算法對(duì)光譜圖像進(jìn)行分塊處理：

其中，x表示的是需要處理的圖像，Ri表示的是在某個(gè)像素點(diǎn)i處從原始圖像x中提取圖像塊（塊的尺寸為×）的矩陣提取操作，xi表示的是將提取出來(lái)的塊重新排列成一個(gè)列向量（尺寸為n），假設(shè)圖像的尺寸為M×N，則總共可以提取出個(gè)不同的圖像塊。

對(duì)圖像進(jìn)行分塊以列向量的形式堆疊成矩陣X，然后對(duì)這個(gè)圖像塊矩陣進(jìn)行排序，即尋找一個(gè)變換矩陣P使得Xp=PX。通常假定對(duì)于干凈圖像Xp是平滑的，Ram I等[17]基于這一假設(shè)提出了分塊排序的降噪算法，用總變分作為“平滑性”的測(cè)量標(biāo)準(zhǔn)：

最小化‖Xp‖TV需要找到一個(gè)最短的路徑使得每個(gè)點(diǎn)只被訪問(wèn)一次，Ram I采用一種近似的解決方案，具體為：從任意點(diǎn)開(kāi)始，每一點(diǎn)xj0到下一點(diǎn)xj1的概率為：

到第二近的點(diǎn)xj2的可能性為：

這里ε是預(yù)先設(shè)定的參數(shù)，α被選取適當(dāng)?shù)闹凳沟胮1+p2=1，xj1和xj2取自未被訪問(wèn)的點(diǎn)集。每一點(diǎn)的搜索范圍限定在以該點(diǎn)為中心的一定大小的方塊中，如果該方塊中沒(méi)有訪問(wèn)的點(diǎn)則在整個(gè)圖像的未訪問(wèn)點(diǎn)中尋找最近鄰點(diǎn)。本文中我們引入這種排序方法進(jìn)行分塊排序。排序完成后用每一點(diǎn)兩側(cè)的小塊拼接成2N×1的向量yi代替原來(lái)的N×1的小塊xi，供下一步重采樣。

2.2 重采樣方法

對(duì)提取出的圖像塊向量進(jìn)行重采樣處理，以獲得更多的樣本數(shù)據(jù)。重采樣的方法可以用以下公式表示：

2.3 權(quán)重計(jì)算

我們可以利用NLM算法中權(quán)重的思想來(lái)對(duì)光譜圖像進(jìn)行降噪處理。NLM[17]算法中像素間的權(quán)重是通過(guò)式(7)來(lái)進(jìn)行計(jì)算的：

其中，h是濾波器參數(shù)，其控制的是濾波器的延滯操作。當(dāng)h→0，降噪后的圖像更加趨近于輸入的噪聲圖像；當(dāng)h→∞，則輸出圖像變得更為模糊。Wi是一個(gè)正則化因子，是所有權(quán)重的總和Wi=∑j∈Iωij。zηi和zηj分別代表的是把塊ηi和ηj展開(kāi)為列向量的強(qiáng)度值，‖zηi-zηj表示的是利用像素間的歐式距離來(lái)進(jìn)行權(quán)重計(jì)算。

在本算法中，我們利用隨機(jī)距離[18]的思想來(lái)選取近似圖像塊。假設(shè)兩個(gè)像素點(diǎn)X和Y分別是服從參數(shù)為λX和λY的泊松隨機(jī)變量，則這兩個(gè)像素點(diǎn)X和Y之間的距離可以表示為

其中，disXY表示的是像素點(diǎn)X和Y之間的距離。

因此，參照式（8），對(duì)于上述通過(guò)重采樣得到的圖像塊樣本矩陣，我們利用隨機(jī)距離來(lái)計(jì)算重采樣樣本矩陣每一列向量數(shù)據(jù)與原始被采集的子塊向量數(shù)據(jù)之間的權(quán)重關(guān)系，即

其中，dis(a, b)用式（8）來(lái)表示，Wi表示的是所有權(quán)重的總和Wi=∑j∈Iωij，λ(ηi)表示的是原始被采集的向量的期望值，λ(Yjxi)表示的是重采樣得到的樣本矩陣Yjxi的每一個(gè)列向量的期望值，γ為指數(shù)的控制因子，設(shè)為常數(shù)。然后對(duì)這些求到的權(quán)重?cái)?shù)據(jù)從大到小地進(jìn)行排序，獲取前M(這里為5)個(gè)最大權(quán)重所對(duì)應(yīng)的樣本矩陣中的列向量，這些列向量就可以被選取為由塊xi所提取的近似圖像塊Dxi，即：

其中，Dxi為選取的近似圖像塊矩陣，為前M個(gè)較大的權(quán)重，Yjxi為塊xi通過(guò)重采樣得到的樣本矩陣。因此，圖像中每個(gè)塊的近似塊經(jīng)過(guò)這樣的選取，可以構(gòu)成一個(gè)最終的近似圖像塊矩陣D，即

2.4 基于指數(shù)分布族的PCA算法

我們假定重采樣得到的相似圖像塊仍然是服從泊松分布的。對(duì)于指數(shù)分布的泊松噪聲圖像數(shù)據(jù)而言，基于NLPCA算法[13]的思想，該噪聲數(shù)據(jù)通?？梢允褂靡韵鹿絹?lái)進(jìn)行分解：

因此，我們的目的是要最小化式（14）來(lái)得到最終的系數(shù)矩陣U*和數(shù)據(jù)主要成分矩陣V*：

因此，最后得到的降噪圖像塊矩陣為：

然后我們?cè)倮檬?3)的逆操作對(duì)求得的圖像塊矩陣進(jìn)行還原，得到最終的降噪圖像Y。為了達(dá)到較好的降噪效果，我們對(duì)第一次降噪后的圖像結(jié)果再次進(jìn)行相同的處理，經(jīng)過(guò)多次迭代，最終達(dá)到理想的效果。

3 實(shí)驗(yàn)結(jié)果

由于我們的算法是針對(duì)小尺寸圖像進(jìn)行處理的，因此，我們256×256的標(biāo)準(zhǔn)Saturn灰度圖像中截取出一塊50×50像素的圖像塊進(jìn)行算法降噪性能的比較。我們之所以考慮小尺寸圖像是由于它們更存在了小樣本和區(qū)塊間圖像特征自相似性不足的問(wèn)題。我們對(duì)這些圖像加入泊松噪聲(其中噪聲峰值分別為peak=0.1以及peak=1)，得到噪聲圖像。然后利用本文方法對(duì)圖像進(jìn)行降噪處理。為了比較我們提出算法的性能優(yōu)劣，實(shí)驗(yàn)中還采用了結(jié)合Anscombe變換[6]的BM3D算法[2]以及NLPCA算法[13]進(jìn)行降噪處理。圖1表示的是不同降噪方法進(jìn)行降噪所得的結(jié)果。圖中(a)到(e)分別表示的是無(wú)噪圖像，噪聲圖像，NLPCA方法的降噪結(jié)果，BM3D算法降噪結(jié)果以及我們算法的降噪結(jié)果。從直觀的視覺(jué)評(píng)判可以看出，我們的算法較其他兩種降噪算法具有更好的降噪性能。

圖1 從標(biāo)準(zhǔn)Saturn圖像中截取出50×50的圖像的降噪處理結(jié)果Fig.1 The denoised results for 50×50 seized by standard Saturn image

為了能夠客觀地衡量算法的降噪性能，表1給出了用峰值信噪比（Peak Signal-to-noise Ratios，PSNR）衡量的性能指標(biāo)。從表1中看出，本方法綜合地優(yōu)于其他兩種算法，特別是在對(duì)小尺寸圖像進(jìn)行降噪處理，更能體現(xiàn)本文算法的優(yōu)勢(shì)。這是因?yàn)槲覀兊乃惴ㄔ谥夭蓸硬僮饕院?，更能充分捕捉圖像本身的統(tǒng)計(jì)特性和所表征的圖像結(jié)果特征。

表1 不同算法降噪結(jié)果的PSNR比較（dB）Table 1 The PSNR (dB) results of different algorithms

此外，我們使用真實(shí)的老鼠圖像來(lái)對(duì)我們的算法性能進(jìn)行評(píng)估。如圖2所示。從圖2中我們可以看出在針對(duì)小尺寸圖像的情況下，相較于其他兩種算法，我們可以得到一個(gè)比較不錯(cuò)的降噪性能。

圖2 截取出的老鼠圖像降噪結(jié)果Fig.2 The denoised results of seized image of mouse

4 總結(jié)

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

[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

國(guó)家自然科學(xué)基金委面上項(xiàng)目（61271320，60872102）；上海交通大學(xué)醫(yī)工交叉基金面上項(xiàng)目(YG2014MS29)

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