沈 帆,李翰林,孫 斌,喻春雨*
(1.南京郵電大學(xué) 電子與光學(xué)工程學(xué)院、微電子學(xué)院,江蘇 南京 210023;2.南京郵電大學(xué) 自動(dòng)化學(xué)院,江蘇 南京 210023)
自1895年X射線被發(fā)現(xiàn)以來(lái),X射線成像技術(shù)被廣泛應(yīng)用于工業(yè)、醫(yī)療等領(lǐng)域[1]。由于泊松噪聲是影響X射線圖像質(zhì)量的主要因素,因此研究泊松噪聲的消除方法及性能提高一直是恢復(fù)X射線圖像質(zhì)量的有效途徑[2-3]。在此方面較為典型的研究有:Selesnick等人提出非高斯雙變量分布,通過(guò)貝葉斯理論導(dǎo)出非線性閾值函數(shù)[4],并結(jié)合隱馬爾科夫模型實(shí)現(xiàn)對(duì)泊松噪聲進(jìn)行濾除[5];Ma,LY等人將圖像的稀疏表示用于圖像恢復(fù),基于像素總變差正則化和數(shù)據(jù)保真度獲得泊松噪聲統(tǒng)計(jì),從而濾除泊松噪聲[6]。此外,近年來(lái)還有對(duì)X射線圖像進(jìn)行Anscombe變換,將圖像中泊松噪聲變換為高斯噪聲,然后利用相對(duì)成熟的高斯噪聲降噪算法對(duì)其進(jìn)行降噪[7],如:Manduca等人提出將Anscombe變換結(jié)合雙邊濾波方法應(yīng)用于低劑量CT投影數(shù)據(jù)降噪[8];畢一鳴等人提出通過(guò)Anscombe變換和三維塊匹配(Block-Matching 3D,BM3D)濾波相結(jié)合進(jìn)行降噪,最后用濾波反投影算法重建圖像[9]。
由于圖像中信號(hào)和噪聲被視為不相關(guān)分量,則基于盲源分離(Blind Source Separation,BSS)理論,將含噪聲圖像看作信號(hào)分量和噪聲分量的組合;而且通過(guò)Anscombe變換可以將泊松噪聲轉(zhuǎn)化為高斯噪聲進(jìn)行消除,由此提出基于Anscombe變換的X射線圖像序列盲源分離降噪方法(以下稱本文降噪算法),其算法思想是:首先取一序列X射線圖像,并利用Anscombe變換將圖像中泊松噪聲變換為高斯噪聲;然后將含噪聲圖像視為信號(hào)分量與噪聲分量的混合,通過(guò)非線性主分量分析(Nonlinear Principal Component Analysis,NLPCA)[10]對(duì)圖像序列進(jìn)行BSS運(yùn)算;在分離分量中取方差較大的進(jìn)行Anscombe逆變換,得到降噪圖像。
首先獲取一個(gè)X射線圖像序列,該序列包含至少2張圖像,然后對(duì)每張圖像進(jìn)行Anscombe變換;將圖像中泊松分布噪聲轉(zhuǎn)化為高斯噪聲。Anscombe變換是Anscombe于1948年提出的一種從泊松分布轉(zhuǎn)到近似高斯分布的非線性變換方法[11-12],如式(1)所示:
(1)
其中:x服從泊松分布,x′是近似服從方差為1的高斯分布。
通過(guò)式(1)產(chǎn)生一組新的含噪聲圖像序列。將每張圖像視為一個(gè)觀測(cè)信號(hào),則該含噪聲圖像序列視為多個(gè)觀測(cè)信號(hào);每個(gè)觀測(cè)信號(hào)是噪聲分量與圖像信號(hào)分量的混合,因此可以用BSS處理將圖像信號(hào)分量從噪聲中分離出來(lái)[13]。分離之后得到的是初步降噪圖像,要通過(guò)Anscombe逆變換才可以得到最終降噪圖像。這里采用NLPCA進(jìn)行BSS處理,是因?yàn)橹鞣至糠治?Principal Component Analysis,PCA)用一系列直線來(lái)描述數(shù)據(jù)中的主要變化趨勢(shì),而NLPCA使用若干開放曲線或閉合曲線來(lái)描述數(shù)據(jù)變化趨勢(shì),更適合描述圖像中信號(hào)和噪聲之間的關(guān)系[14],其優(yōu)化條件如式(2)所示:
(2)
其中:wiTx′為非線性因子,g(x′)是非線性函數(shù),其表達(dá)如式(3)所示:
g(x′)=x-G(x′)+G(ν),
(3)
其中:ν為高斯隨機(jī)向量,取G(x)=tanh(ax),a為一常數(shù)。
假設(shè)W是一個(gè)正交矩陣,令S=Wx′,則用于優(yōu)化的指標(biāo)函數(shù)如式(4)所示:
(4)
其中:算法采用非線性遞歸最小二乘學(xué)習(xí)規(guī)則[15],它根據(jù)輸入精度數(shù)據(jù)自動(dòng)決定適合的學(xué)習(xí)速率參數(shù),具有比隨機(jī)梯度方法收斂快的優(yōu)勢(shì)[16]。
圖1 本文降噪算法流程圖
圖1為本文降噪算法的流程圖,算法流程依次為:對(duì)X射線圖像序列中圖像進(jìn)行Anscombe變換;通過(guò)NLPCA對(duì)Anscombe變換后含噪聲圖像進(jìn)行BSS處理,得到初步降噪圖像;進(jìn)行Anscombe逆變換,得到最終降噪圖像。由此制定算法步驟如下:
Step 1:利用X射線成像設(shè)備拍攝一組圖像序列,所拍攝目標(biāo)與場(chǎng)景相對(duì)靜止,采樣張數(shù)為n;
Step 2:利用Anscombe變換將X射線圖像中泊松噪聲轉(zhuǎn)化為高斯噪聲;
Step 3:將經(jīng)過(guò)Anscombe變換后的二維圖像轉(zhuǎn)換為一維數(shù)組,得到n個(gè)觀測(cè)變量;
Step 4:對(duì)觀測(cè)到的混合信號(hào)進(jìn)行中心化處理,即去均值處理[17];
Step 5:取t=0,生成隨機(jī)初始矩陣為W(0);
Step 6:按照非線性遞歸最小二乘規(guī)則迭代;
Step 7:疊加次數(shù)t=t+1,判斷此次是否滿足|J(W(t))-J(W(t-1))| Step 8:將分離得到的圖像信號(hào)再通過(guò)Anscombe逆變換,得到最終降噪圖像。 本文分別通過(guò)Shepp-Logan頭模型圖像和自制X射線機(jī)拍攝的圖像對(duì)本文降噪算法進(jìn)行性能分析。采用Shepp-Logan頭模型圖像是因?yàn)橛伤梢噪S機(jī)生成任意張數(shù)、任意噪聲強(qiáng)度的含噪聲圖像,這有利于對(duì)本文降噪算法的降噪局限進(jìn)行分析;采用真實(shí)低質(zhì)量X射線圖像來(lái)分析降噪方法可以驗(yàn)證算法的實(shí)用性。本文降噪算法的計(jì)算平臺(tái)為Intel(R) Core(TM)i7-6700 CPU、@3.40 GHz 3.40 GHz四核處理器、8 GB內(nèi)存的PC機(jī)。 如圖2(a),Sheep-Logan頭模型被用來(lái)表示一個(gè)腦部斷層圖像,它由10個(gè)位置大小、方向、密度各不相同的橢圓組成,通過(guò)不同的橢圓來(lái)表征不同形狀,利用不同灰度值模擬不同組織的衰減系數(shù)[18]。 圖2 原始圖、含噪聲圖及降噪圖 圖2(a)為Sheep-Logan原始圖像;圖2(b)為含泊松噪聲圖像,其噪聲是以圖2(a)中像素值為均值隨機(jī)生成的,其峰值信噪比(Peak Signal to Noise Ratio, PSNR)為28.289 4 dB、結(jié)構(gòu)相似性(Structural Similarity Index, SSIM)為0.700 7;圖2(c)~圖2(i)對(duì)應(yīng)采樣序列中含不同圖像數(shù)目時(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時(shí),降噪圖像的PSNR和SSIM值明顯提高,PSNR由31.859 8 dB提高到37.267 8 dB、SSIM由0.790 2提高到0.963 8。當(dāng)采樣序列中n數(shù)值增加到50時(shí),本文降噪算法通過(guò)13次迭代即可以完成收斂,運(yùn)行時(shí)間為0.738 4 s;繼續(xù)增加采樣張數(shù),算法運(yùn)行時(shí)間增加,降噪效果無(wú)明顯改善。 圖3 降噪圖像PSNR、SSIM隨序列中圖像張數(shù)的變化 相應(yīng)于表1的代表性數(shù)據(jù),圖3截取本文降噪算法對(duì)Sheep-Logan頭模型的降噪效果隨采樣序列中圖像張數(shù)的變化規(guī)律曲線中[2,120]段:在圖像采樣數(shù)目處于[2, 20]時(shí),隨著圖像數(shù)目增加,降噪效果提升明顯;圖像采樣數(shù)目處于[20, 50],圖像的PSNR值以及SSIM值提升緩慢,圖像采樣數(shù)目為50,降噪效果最佳;圖像采樣數(shù)目處于[50, +∞),降噪效果無(wú)明顯提升,但可以通過(guò)改變算法精度完善。 圖4 不同降噪算法的降噪效果圖 此外,圖4給出本文降噪算法同幾種常見(jiàn)消除泊松噪聲的算法降噪效果比較。圖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。在視覺(jué)效果上比較:圖4(b)噪聲殘留較多;圖4(c)細(xì)節(jié)輪廓不全;圖4(d)圖像亮度較差;圖4(e)和圖4(f)在圖像亮度和細(xì)節(jié)輪廓保留上有優(yōu)勢(shì),且圖4(f)效果最優(yōu)。 表2 不同降噪算法性能對(duì)比 表2為對(duì)應(yīng)圖4的算法性能參數(shù)對(duì)比,對(duì)比評(píng)價(jià)參數(shù):本文降噪算法在降噪效果上最優(yōu),且在運(yùn)行時(shí)間上也具有明顯優(yōu)勢(shì)。 在使用真實(shí)X射線圖像對(duì)本文降噪算法進(jìn)行性能分析時(shí),選用實(shí)驗(yàn)組自制的X射線機(jī)拍攝9幀像素為768×576的X射線圖像組成圖像序列。由于無(wú)“干凈”原圖像做降噪對(duì)比,這里采用信噪比(Signal Noise Ratio,SNR)、信息熵(Entropy)、標(biāo)準(zhǔn)差(Standard Deviation, SD)和運(yùn)行時(shí)間(Runtime)作評(píng)價(jià)參數(shù),與參考文獻(xiàn)[19]的SVD+Sequence降噪方法做對(duì)比。 圖5 兩種算法在不同采樣圖像數(shù)目下的降噪效果對(duì)比 圖5中,圖5(a)為采用本文降噪算法的降噪效果及其邊緣檢測(cè)圖,圖5(b)為采用SVD+Sequence降噪算法的降噪效果及相應(yīng)邊緣檢測(cè)圖。 圖6 SNR,SD,Entropy和Runtime上的對(duì)比 圖6給出了本文降噪算法和參考文獻(xiàn)[19]降噪算法的對(duì)比曲線,橫坐標(biāo)均為序列中采樣圖像數(shù)目。圖6(a)為SNR對(duì)比,由此可見(jiàn)當(dāng)序列中圖像數(shù)目由2增加到9時(shí),本文降噪算法的SNR值由30.355 1 dB提高到65.949 8 dB,高于SVD+Sequence算法的SNR值;圖6(b)為SD對(duì)比,可見(jiàn)本文降噪算法的降噪圖像SD較低,說(shuō)明其受像素值干擾較??;圖6(c)為Entropy對(duì)比,可見(jiàn)當(dāng)序列中圖像數(shù)目由2增加到9時(shí),本文降噪算法圖像Entropy從5.923 5提高到6.022 3,值優(yōu)于SVD+Sequence算法,說(shuō)明本文降噪算法更有效保留細(xì)節(jié)信息;圖6(d)為Runtime對(duì)比,可見(jiàn)本文降噪算法運(yùn)行時(shí)間較長(zhǎng),采用9張X射線圖像降噪時(shí),迭代16次,Runtime為1.311 s,同樣滿足實(shí)時(shí)觀測(cè)需要。 由圖5和圖6得出結(jié)論:當(dāng)參與降噪的圖像張數(shù)增加時(shí),兩種算法的降噪效果均有改善且圖像邊緣更完整;本文降噪算法在降噪效果及邊緣細(xì)節(jié)保留上更優(yōu);本文降噪算法在運(yùn)算時(shí)間上不占優(yōu)勢(shì)。 圖7 其他常用降噪算法的降噪效果圖 圖7給出本文降噪算法同其他幾種常用降噪算法的降噪效果對(duì)比,這里采樣序列中圖像數(shù)目為9。圖7(a),7(d)分別為Wavelet降噪圖及對(duì)應(yīng)邊緣檢測(cè)圖;圖7(b),7(e)為Bilateral filter降噪圖及對(duì)應(yīng)邊緣檢測(cè)圖;圖7(c),7(f)為Anscombe+BM3D降噪圖及對(duì)應(yīng)邊緣檢測(cè)圖。由對(duì)比看出,相比較圖5(a)中本文降噪算法(n=9),常見(jiàn)的這3種算法得到的降噪圖像中噪聲殘留更多,圖像細(xì)節(jié)、邊緣保留較差。這說(shuō)明:當(dāng)序列中取適當(dāng)數(shù)量含噪聲圖像時(shí),本文降噪算法更具性能優(yōu)勢(shì),不僅可以有效分離噪聲,還可以有效保留圖像的輪廓邊緣和細(xì)節(jié)信息。 表3為圖7所示圖像的SNR,Entropy,SD以及Runtime。圖像SNR值高說(shuō)明算法對(duì)于噪聲濾除更有效;Entropy高則說(shuō)明降噪后圖像中細(xì)節(jié)信息保存較多;圖像SD大說(shuō)明像素值受干擾較大。分析數(shù)據(jù)說(shuō)明:本文降噪算法在降噪和信息保留上均優(yōu)于其他算法,同時(shí)算法運(yùn)行時(shí)間也滿足實(shí)時(shí)需要。 表3 常用降噪算法降噪性能對(duì)比 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 本研究介紹一種基于Anscombe變換的X射線圖像序列盲源分離降噪方法,并分別利用Sheep-Logan模型圖像和真實(shí)X射線圖像對(duì)該降噪方法的性能進(jìn)行評(píng)價(jià)分析。分析結(jié)果表明:隨著圖像序列中采樣張數(shù)增加,本文算法降噪效果顯著提高,并在采樣張數(shù)達(dá)到50的時(shí)候,降噪性能達(dá)到最優(yōu);增加圖像序列中采樣張數(shù)會(huì)提高降噪性能,但也會(huì)增加算法運(yùn)行時(shí)間,需要根據(jù)實(shí)際情況找到運(yùn)行時(shí)間和降噪性能的平衡;本文降噪算法性能可以通過(guò)GPU得以大大改善,是一種實(shí)用的降噪方法。3 本文降噪算法性能分析
3.1 采用Sheep-Logan頭模圖像對(duì)本文降噪算法進(jìn)行性能分析
3.2 真實(shí)X射線圖像對(duì)本文降噪算法的性能分析
4 結(jié) 論