沈 帆，李翰林，孫 斌，喻春雨*

(1.南京郵電大學(xué) 電子與光學(xué)工程學(xué)院、微電子學(xué)院，江蘇 南京 210023；2.南京郵電大學(xué) 自動化學(xué)院，江蘇 南京 210023)

1 引 言

自1895年X射線被發(fā)現(xiàn)以來，X射線成像技術(shù)被廣泛應(yīng)用于工業(yè)、醫(yī)療等領(lǐng)域[1]。由于泊松噪聲是影響X射線圖像質(zhì)量的主要因素，因此研究泊松噪聲的消除方法及性能提高一直是恢復(fù)X射線圖像質(zhì)量的有效途徑[2-3]。在此方面較為典型的研究有：Selesnick等人提出非高斯雙變量分布，通過貝葉斯理論導(dǎo)出非線性閾值函數(shù)[4]，并結(jié)合隱馬爾科夫模型實現(xiàn)對泊松噪聲進(jìn)行濾除[5]；Ma，LY等人將圖像的稀疏表示用于圖像恢復(fù)，基于像素總變差正則化和數(shù)據(jù)保真度獲得泊松噪聲統(tǒng)計，從而濾除泊松噪聲[6]。此外，近年來還有對X射線圖像進(jìn)行Anscombe變換，將圖像中泊松噪聲變換為高斯噪聲，然后利用相對成熟的高斯噪聲降噪算法對其進(jìn)行降噪[7]，如：Manduca等人提出將Anscombe變換結(jié)合雙邊濾波方法應(yīng)用于低劑量CT投影數(shù)據(jù)降噪[8]；畢一鳴等人提出通過Anscombe變換和三維塊匹配(Block-Matching 3D，BM3D)濾波相結(jié)合進(jìn)行降噪，最后用濾波反投影算法重建圖像[9]。

由于圖像中信號和噪聲被視為不相關(guān)分量，則基于盲源分離(Blind Source Separation，BSS)理論，將含噪聲圖像看作信號分量和噪聲分量的組合；而且通過Anscombe變換可以將泊松噪聲轉(zhuǎn)化為高斯噪聲進(jìn)行消除，由此提出基于Anscombe變換的X射線圖像序列盲源分離降噪方法(以下稱本文降噪算法)，其算法思想是：首先取一序列X射線圖像，并利用Anscombe變換將圖像中泊松噪聲變換為高斯噪聲；然后將含噪聲圖像視為信號分量與噪聲分量的混合，通過非線性主分量分析(Nonlinear Principal Component Analysis，NLPCA)[10]對圖像序列進(jìn)行BSS運算；在分離分量中取方差較大的進(jìn)行Anscombe逆變換，得到降噪圖像。

2 本文降噪算法介紹

首先獲取一個X射線圖像序列，該序列包含至少2張圖像，然后對每張圖像進(jìn)行Anscombe變換；將圖像中泊松分布噪聲轉(zhuǎn)化為高斯噪聲。Anscombe變換是Anscombe于1948年提出的一種從泊松分布轉(zhuǎn)到近似高斯分布的非線性變換方法[11-12]，如式(1)所示：

(1)

其中：x服從泊松分布，x′是近似服從方差為1的高斯分布。

通過式(1)產(chǎn)生一組新的含噪聲圖像序列。將每張圖像視為一個觀測信號，則該含噪聲圖像序列視為多個觀測信號；每個觀測信號是噪聲分量與圖像信號分量的混合，因此可以用BSS處理將圖像信號分量從噪聲中分離出來[13]。分離之后得到的是初步降噪圖像，要通過Anscombe逆變換才可以得到最終降噪圖像。這里采用NLPCA進(jìn)行BSS處理，是因為主分量分析(Principal Component Analysis，PCA)用一系列直線來描述數(shù)據(jù)中的主要變化趨勢，而NLPCA使用若干開放曲線或閉合曲線來描述數(shù)據(jù)變化趨勢，更適合描述圖像中信號和噪聲之間的關(guān)系[14]，其優(yōu)化條件如式(2)所示：

(2)

其中：wiTx′為非線性因子，g(x′)是非線性函數(shù)，其表達(dá)如式(3)所示：

g(x′)=x-G(x′)+G(ν),

(3)

其中：ν為高斯隨機向量，取G(x)=tanh(ax)，a為一常數(shù)。

假設(shè)W是一個正交矩陣，令S=Wx′，則用于優(yōu)化的指標(biāo)函數(shù)如式(4)所示：

(4)

其中：算法采用非線性遞歸最小二乘學(xué)習(xí)規(guī)則[15]，它根據(jù)輸入精度數(shù)據(jù)自動決定適合的學(xué)習(xí)速率參數(shù)，具有比隨機梯度方法收斂快的優(yōu)勢[16]。

圖1 本文降噪算法流程圖

圖1為本文降噪算法的流程圖，算法流程依次為：對X射線圖像序列中圖像進(jìn)行Anscombe變換；通過NLPCA對Anscombe變換后含噪聲圖像進(jìn)行BSS處理，得到初步降噪圖像；進(jìn)行Anscombe逆變換，得到最終降噪圖像。由此制定算法步驟如下：

Step 1：利用X射線成像設(shè)備拍攝一組圖像序列，所拍攝目標(biāo)與場景相對靜止，采樣張數(shù)為n；

Step 2：利用Anscombe變換將X射線圖像中泊松噪聲轉(zhuǎn)化為高斯噪聲；

Step 3：將經(jīng)過Anscombe變換后的二維圖像轉(zhuǎn)換為一維數(shù)組，得到n個觀測變量；

Step 4：對觀測到的混合信號進(jìn)行中心化處理，即去均值處理[17]；

Step 5：取t=0，生成隨機初始矩陣為W(0)；

Step 6：按照非線性遞歸最小二乘規(guī)則迭代；

Step 7：疊加次數(shù)t=t+1，判斷此次是否滿足|J(W(t))-J(W(t-1))|

Step 8：將分離得到的圖像信號再通過Anscombe逆變換，得到最終降噪圖像。

3 本文降噪算法性能分析

本文分別通過Shepp-Logan頭模型圖像和自制X射線機拍攝的圖像對本文降噪算法進(jìn)行性能分析。采用Shepp-Logan頭模型圖像是因為由它可以隨機生成任意張數(shù)、任意噪聲強度的含噪聲圖像，這有利于對本文降噪算法的降噪局限進(jìn)行分析；采用真實低質(zhì)量X射線圖像來分析降噪方法可以驗證算法的實用性。本文降噪算法的計算平臺為Intel(R) Core(TM)i7-6700 CPU、@3.40 GHz 3.40 GHz四核處理器、8 GB內(nèi)存的PC機。

3.1 采用Sheep-Logan頭模圖像對本文降噪算法進(jìn)行性能分析

如圖2(a)，Sheep-Logan頭模型被用來表示一個腦部斷層圖像，它由10個位置大小、方向、密度各不相同的橢圓組成，通過不同的橢圓來表征不同形狀，利用不同灰度值模擬不同組織的衰減系數(shù)[18]。

圖2 原始圖、含噪聲圖及降噪圖

圖2(a)為Sheep-Logan原始圖像；圖2(b)為含泊松噪聲圖像，其噪聲是以圖2(a)中像素值為均值隨機生成的，其峰值信噪比(Peak Signal to Noise Ratio, PSNR)為28.289 4 dB、結(jié)構(gòu)相似性(Structural Similarity Index, SSIM)為0.700 7；圖2(c)～圖2(i)對應(yīng)采樣序列中含不同圖像數(shù)目時的分離降噪圖像；n表示采樣序列中圖像數(shù)目。由圖2看出：隨著采樣序列中圖像數(shù)目的增加，降噪圖像的質(zhì)量得到明顯改善。

表1 降噪圖像性能隨采樣序列中圖像張數(shù)的變化

Tab.1 Denoising performance varies with number of images in sequence

采樣數(shù)/張PSNR/dBSSIMRuntime/s231.859 80.790 20.122 91036.105 30.932 70.279 12036.795 60.951 90.399 63037.045 70.958 20.551 25037.267 80.963 80.738 48037.184 50.962 41.371 412036.829 40.958 61.943 6

表1所示為降噪圖像的PSNR、SSIM和Runtime。當(dāng)采樣序列中n數(shù)值從2增加到50時，降噪圖像的PSNR和SSIM值明顯提高，PSNR由31.859 8 dB提高到37.267 8 dB、SSIM由0.790 2提高到0.963 8。當(dāng)采樣序列中n數(shù)值增加到50時，本文降噪算法通過13次迭代即可以完成收斂，運行時間為0.738 4 s；繼續(xù)增加采樣張數(shù)，算法運行時間增加，降噪效果無明顯改善。

圖3 降噪圖像PSNR、SSIM隨序列中圖像張數(shù)的變化

相應(yīng)于表1的代表性數(shù)據(jù)，圖3截取本文降噪算法對Sheep-Logan頭模型的降噪效果隨采樣序列中圖像張數(shù)的變化規(guī)律曲線中[2,120]段：在圖像采樣數(shù)目處于[2, 20]時，隨著圖像數(shù)目增加，降噪效果提升明顯；圖像采樣數(shù)目處于[20, 50]，圖像的PSNR值以及SSIM值提升緩慢，圖像采樣數(shù)目為50，降噪效果最佳；圖像采樣數(shù)目處于[50, +∞)，降噪效果無明顯提升，但可以通過改變算法精度完善。

圖4 不同降噪算法的降噪效果圖

此外，圖4給出本文降噪算法同幾種常見消除泊松噪聲的算法降噪效果比較。圖4(a)為含噪聲圖像；圖4(b)為軟閾值小波(Wavelet)降噪圖像；圖4(c)為雙邊濾波(Bilateral filter)降噪圖像；圖4(d)為Anscombe+BM3D降噪圖像；圖4(e)為參考文獻(xiàn)[19]的SVD+Sequence降噪圖像[19]，序列中圖像數(shù)目為9；圖4(f)為本文降噪算法的降噪圖像，序列中圖像數(shù)目同樣取9。在視覺效果上比較：圖4(b)噪聲殘留較多；圖4(c)細(xì)節(jié)輪廓不全；圖4(d)圖像亮度較差；圖4(e)和圖4(f)在圖像亮度和細(xì)節(jié)輪廓保留上有優(yōu)勢，且圖4(f)效果最優(yōu)。

表2 不同降噪算法性能對比

表2為對應(yīng)圖4的算法性能參數(shù)對比，對比評價參數(shù)：本文降噪算法在降噪效果上最優(yōu)，且在運行時間上也具有明顯優(yōu)勢。

3.2 真實X射線圖像對本文降噪算法的性能分析

在使用真實X射線圖像對本文降噪算法進(jìn)行性能分析時，選用實驗組自制的X射線機拍攝9幀像素為768×576的X射線圖像組成圖像序列。由于無“干凈”原圖像做降噪對比，這里采用信噪比(Signal Noise Ratio，SNR)、信息熵(Entropy)、標(biāo)準(zhǔn)差(Standard Deviation, SD)和運行時間(Runtime)作評價參數(shù)，與參考文獻(xiàn)[19]的SVD+Sequence降噪方法做對比。

圖5 兩種算法在不同采樣圖像數(shù)目下的降噪效果對比

圖5中，圖5(a)為采用本文降噪算法的降噪效果及其邊緣檢測圖，圖5(b)為采用SVD+Sequence降噪算法的降噪效果及相應(yīng)邊緣檢測圖。

圖6 SNR，SD，Entropy和Runtime上的對比

圖6給出了本文降噪算法和參考文獻(xiàn)[19]降噪算法的對比曲線，橫坐標(biāo)均為序列中采樣圖像數(shù)目。圖6(a)為SNR對比，由此可見當(dāng)序列中圖像數(shù)目由2增加到9時，本文降噪算法的SNR值由30.355 1 dB提高到65.949 8 dB，高于SVD+Sequence算法的SNR值；圖6(b)為SD對比，可見本文降噪算法的降噪圖像SD較低，說明其受像素值干擾較?。粓D6(c)為Entropy對比，可見當(dāng)序列中圖像數(shù)目由2增加到9時，本文降噪算法圖像Entropy從5.923 5提高到6.022 3，值優(yōu)于SVD+Sequence算法，說明本文降噪算法更有效保留細(xì)節(jié)信息；圖6(d)為Runtime對比，可見本文降噪算法運行時間較長，采用9張X射線圖像降噪時，迭代16次，Runtime為1.311 s，同樣滿足實時觀測需要。

由圖5和圖6得出結(jié)論：當(dāng)參與降噪的圖像張數(shù)增加時，兩種算法的降噪效果均有改善且圖像邊緣更完整；本文降噪算法在降噪效果及邊緣細(xì)節(jié)保留上更優(yōu)；本文降噪算法在運算時間上不占優(yōu)勢。

圖7 其他常用降噪算法的降噪效果圖

圖7給出本文降噪算法同其他幾種常用降噪算法的降噪效果對比，這里采樣序列中圖像數(shù)目為9。圖7(a)，7(d)分別為Wavelet降噪圖及對應(yīng)邊緣檢測圖；圖7(b)，7(e)為Bilateral filter降噪圖及對應(yīng)邊緣檢測圖；圖7(c)，7(f)為Anscombe+BM3D降噪圖及對應(yīng)邊緣檢測圖。由對比看出，相比較圖5(a)中本文降噪算法(n=9)，常見的這3種算法得到的降噪圖像中噪聲殘留更多，圖像細(xì)節(jié)、邊緣保留較差。這說明：當(dāng)序列中取適當(dāng)數(shù)量含噪聲圖像時，本文降噪算法更具性能優(yōu)勢，不僅可以有效分離噪聲，還可以有效保留圖像的輪廓邊緣和細(xì)節(jié)信息。

表3為圖7所示圖像的SNR，Entropy，SD以及Runtime。圖像SNR值高說明算法對于噪聲濾除更有效；Entropy高則說明降噪后圖像中細(xì)節(jié)信息保存較多；圖像SD大說明像素值受干擾較大。分析數(shù)據(jù)說明：本文降噪算法在降噪和信息保留上均優(yōu)于其他算法，同時算法運行時間也滿足實時需要。

表3 常用降噪算法降噪性能對比

Tab.3 Denoising performance comparison of common algorithms

降噪算法SNR/dBEntropySDRuntime/sWavelet55.304 45.959 870.7960.299Bilateral filter47.274 35.903 470.8151.625Anscombe+BM3D51.377 55.962 571.0354.596SVD+Sequence[19]47.233 75.863 670.8340.381Proposed denoising65.94986.022 370.7821.311

4 結(jié) 論

本研究介紹一種基于Anscombe變換的X射線圖像序列盲源分離降噪方法，并分別利用Sheep-Logan模型圖像和真實X射線圖像對該降噪方法的性能進(jìn)行評價分析。分析結(jié)果表明：隨著圖像序列中采樣張數(shù)增加，本文算法降噪效果顯著提高，并在采樣張數(shù)達(dá)到50的時候，降噪性能達(dá)到最優(yōu)；增加圖像序列中采樣張數(shù)會提高降噪性能，但也會增加算法運行時間，需要根據(jù)實際情況找到運行時間和降噪性能的平衡；本文降噪算法性能可以通過GPU得以大大改善，是一種實用的降噪方法。