基于殘差學(xué)習(xí)的CT稀疏重建偽影校正

張鵬程，謝世朋

(南京郵電大學(xué) 通信與信息工程學(xué)院，江蘇 南京 210003)

0 引 言

輻射劑量問題已成為CT領(lǐng)域亟待解決的熱點(diǎn)問題，當(dāng)前的研究主要集中在降低電流和減少投影數(shù)方面。降低電流的方法，雖然降低了X射線的輻射量，但同時(shí)引入了較多的統(tǒng)計(jì)噪聲。為了去除低劑量CT中的統(tǒng)計(jì)噪聲，L.L.Chen等提出了使用改進(jìn)的BM3D算法來提高圖像的質(zhì)量[1]，而SayedMasoud Hashemi等采用了NCRE的迭代方法減小低劑量CT中的噪聲[2]。減少投影數(shù)的方法可以減少CT中的大量輻射，但同時(shí)給圖像帶來了嚴(yán)重的偽影。為了去除重建中的偽影，Shaoting Zhang等利用字典學(xué)習(xí)的方法稀疏重建[3]，張瀚銘等利用非局部的TV正則化稀疏重建[4]，張俊峰等利用伽瑪正則化稀疏重建[5]，在圖像去偽影方面取得了良好的效果。由于數(shù)據(jù)噪聲的存在和不理想的先驗(yàn)懲罰項(xiàng)的選取，使得基于TV正則化的重建結(jié)果存在梯形偽影和片狀偽影。在臨床診斷過程中，片狀偽影會(huì)影響到診斷的誤判和正確性。其它稀疏重建的方法都需要大量的正則化約束，且恢復(fù)出來的圖像有較低的PSNR(峰值信噪比)和SSIM(結(jié)構(gòu)相似度)值，仍然存在較嚴(yán)重的偽影。卷積神經(jīng)網(wǎng)絡(luò)(CNN)已經(jīng)在圖像去噪、圖像識(shí)別和圖像分類領(lǐng)域取得了很好的效果且誤差率都很小，在醫(yī)學(xué)圖像處理中得到了廣泛地應(yīng)用。針對(duì)這一現(xiàn)象，本文提出了基于殘差學(xué)習(xí)的CT稀疏重建偽影校正的方法，能更好估計(jì)這些偽影的特征達(dá)到去偽影的目的。

1 方 法

本文提出了基于殘差學(xué)習(xí)的CT稀疏重建偽影校正方法，首先通過FBP(濾波反投影)方法對(duì)CT產(chǎn)生的稀疏投影數(shù)據(jù)進(jìn)行重建，重建后的圖像帶有嚴(yán)重的偽影；然后通過建立殘差神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)對(duì)偽影的特征進(jìn)行學(xué)習(xí)，得到偽影圖；最后通過稀疏重建的圖與偽影的圖做殘差，恢復(fù)出清晰的CT圖像。

1.1 FBP重建

濾波反投影重建(filtered back-projection，F(xiàn)BP)算法具有重建速度快、易于硬件實(shí)現(xiàn)等優(yōu)點(diǎn),當(dāng)投影數(shù)據(jù)比較完備時(shí),能快速重建出高質(zhì)量圖像,因此它是目前應(yīng)用最為廣泛的一種重建算法。對(duì)于平行束CT系統(tǒng)，遵循bear定律。這個(gè)投影的過程實(shí)質(zhì)上是一個(gè)拉東變化的過程。探測(cè)器的距離為t，投影角度為θ，公式可以表示為

(1)

f(x,y)代表當(dāng)前的圖像，當(dāng)t=xcosθ+ysinθ，采樣投影數(shù)據(jù)充足那么FBP則有

(2)

其中，w表示斜坡濾波器，Pθ(w)表示沿著探測(cè)器方向的一維傅里葉變換，wPθ(w)就是濾波的投影。第二重積分本質(zhì)上是‘濾波反投影’，即對(duì)濾波投影的數(shù)據(jù)進(jìn)行傅里葉反變換，然后累積得到反投影的圖像。

使用CT對(duì)患者進(jìn)行不同投影數(shù)的掃描，產(chǎn)生60個(gè)(3600圓軌道掃描，掃描間隔為60)、90個(gè)(3600圓軌道掃描，掃描間隔為40)、120個(gè)投影(3600圓軌道掃描，掃描間隔為30)和全投影(3600圓軌道掃描，掃描間隔為10)的數(shù)據(jù)，而后對(duì)這些數(shù)據(jù)進(jìn)行FBP反投影重建。由CT掃描患者肺部切片全投影成像和稀疏重建成像如圖1所示。

圖1 全投影和稀疏投影重建

FBP稀疏重建的CT圖像存在明顯的條狀偽影，投影數(shù)減少能夠降低X射線對(duì)人體造成的輻射危害，但同時(shí)帶來了偽影的增加，為了去除偽影引入了基于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)。

1.2 基于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)

由于深度學(xué)習(xí)在學(xué)習(xí)圖像噪聲特點(diǎn)方面取得了很好的效果，張凱等[6]利用殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)去高斯噪聲、解決SISR和JPEG圖像塊效應(yīng)等方面較傳統(tǒng)的方法效果更好。另外，Yo Seob Han等[7]根據(jù)張凱的思路提出了基于殘差學(xué)習(xí)的壓縮感知CT重建方法，該方法在去稀疏重建偽影方面同樣取得了很好的效果。基于二者的思想，本文利用基于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)對(duì)稀疏重建的圖像進(jìn)行學(xué)習(xí)，然后通過稀疏重建的圖像減去學(xué)習(xí)得到的偽影圖像，最后恢復(fù)出清晰的校正圖像。殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)的輸入為稀疏重建圖像f=x+n。在深度學(xué)習(xí)卷積網(wǎng)絡(luò)當(dāng)中，典型的去噪模型如MLP[8]和CSF[9]旨在學(xué)習(xí)一個(gè)映射函數(shù)F(f)=x預(yù)測(cè)一個(gè)清晰的圖像。對(duì)于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)，訓(xùn)練的是一個(gè)映射μ(f)≈n函數(shù)。本文提出的基于殘差的卷積神經(jīng)網(wǎng)絡(luò)旨在解決以下問題[8]

(3)

(4)

(5)

該殘差學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)即為本文提出的去偽影校正方法，不同于何愷明等[10]提出的用多個(gè)殘差單元預(yù)測(cè)殘差圖像，本文只使用單個(gè)殘差單元預(yù)測(cè)殘差圖像。由殘差學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)得到的殘差圖即為偽影特征圖，然后根據(jù)x=f-n，得到最終清晰的圖像，實(shí)現(xiàn)過程如圖3所示。

圖2 殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

圖3 稀疏重建圖中去偽影恢復(fù)CT的過程

1.2.1 殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)的具體結(jié)構(gòu)

該網(wǎng)絡(luò)的深度為D，它有3種不同的網(wǎng)絡(luò)層，分別為：卷積網(wǎng)絡(luò)+ReLU[11](修正線性單元)、卷積網(wǎng)絡(luò)+批規(guī)范化(batch normalization，BN)和卷積網(wǎng)絡(luò)層。第一層為卷積網(wǎng)絡(luò)和ReLU構(gòu)成的網(wǎng)絡(luò)層，它由64個(gè)3*3*1的濾波器組成，用來產(chǎn)生64個(gè)特征映射。這一層加入ReLU，目的是在卷積網(wǎng)絡(luò)中引入非線性，因?yàn)樯窠?jīng)網(wǎng)絡(luò)學(xué)習(xí)的實(shí)際數(shù)據(jù)是非線性的，卷積是一個(gè)線性操作——元素級(jí)別的矩陣相乘和相加，因此需要通過使用非線性函數(shù)ReLU來引入非線性。然而經(jīng)實(shí)踐證明，ReLU具備引導(dǎo)適度稀疏的能力，訓(xùn)練后的網(wǎng)絡(luò)完全具備適度的稀疏性，而且訓(xùn)練后的可視化效果和傳統(tǒng)預(yù)處理的效果很相似。正因?yàn)樗羞@樣的作用使得ReLU的速度非常快，而且精確度很高。第二層到D-1層，在卷積網(wǎng)絡(luò)和修正線性單元的基礎(chǔ)上引入了BN(批規(guī)范)[12]，使用64個(gè)3*3*64的濾波器，將BN加在卷積和ReLU之間。深層神經(jīng)網(wǎng)絡(luò)在做非線性變換前激活輸入值隨著網(wǎng)絡(luò)深度的加深，分布逐漸發(fā)生偏移或者變動(dòng)，逐漸往非線性函數(shù)取值區(qū)間的上下限兩端靠近，這導(dǎo)致后向傳播時(shí)低層神經(jīng)網(wǎng)絡(luò)梯度的消失，進(jìn)而使得深層神經(jīng)網(wǎng)絡(luò)的訓(xùn)練速度收斂越來越慢。BN通過一定的規(guī)范化手段，把每層神經(jīng)網(wǎng)絡(luò)中任意神經(jīng)元輸入值的分布強(qiáng)行拉回到均值為0、方差為1的標(biāo)準(zhǔn)正態(tài)分布中(即把越來越偏的分布強(qiáng)制拉回比較標(biāo)準(zhǔn)的分布中)，使得激活輸入值落在非線性函數(shù)對(duì)輸入比較敏感的區(qū)域。這樣輸入的小變化引起損失函數(shù)較大的變化，梯度隨之變大，于是避免了梯度消失的問題，同時(shí)學(xué)習(xí)收斂速度加快，大大提高了訓(xùn)練速度。第三層為一個(gè)卷積層，使用3*3*64的濾波器用來重建輸出結(jié)果。

1.2.2 優(yōu)化加速

在訓(xùn)練中優(yōu)化過程一般采用梯度下降法，梯度下降法有隨機(jī)梯度下降法(stochastic gradient descent，SGD)和Adam算法等。其中Adam算法根據(jù)損失函數(shù)對(duì)每個(gè)參數(shù)梯度的一階矩估計(jì)和二階矩估計(jì)進(jìn)行動(dòng)態(tài)調(diào)整，使其調(diào)整到一個(gè)合適的學(xué)習(xí)速率。Adam算法在每次迭代中參數(shù)的學(xué)習(xí)步長都有一個(gè)確定的范圍，不會(huì)因?yàn)楹艽蟮奶荻葘?dǎo)致學(xué)習(xí)步長的增大，參數(shù)的值也比較穩(wěn)定。批量梯度下降算法的缺點(diǎn)是每次都會(huì)使用全部訓(xùn)練樣本，而且每次都使用完全相同的樣本集，這會(huì)使得計(jì)算相對(duì)的冗雜。相反隨機(jī)梯度下降算法的實(shí)現(xiàn)是可以在線更新的，每次只隨機(jī)選擇一個(gè)樣本來更新模型參數(shù)，學(xué)習(xí)速度是非?？斓摹Ｔ谏窠?jīng)網(wǎng)絡(luò)訓(xùn)練過程當(dāng)中，SGD和Adam優(yōu)化方法的引入使網(wǎng)絡(luò)結(jié)構(gòu)訓(xùn)練得到了更好的結(jié)果，其中Adam算法在圖像高斯去噪過程中效果要強(qiáng)與SGD算法[7]。

2 結(jié) 果

2.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)的數(shù)據(jù)由南京某醫(yī)院提供，來自10個(gè)患者肺部的CT掃描。實(shí)驗(yàn)設(shè)備為臨床使用的西門子SOMATOM Sensation 16排CT，曝光量為100 mAs，管電壓為120 kVp,射線源距離探測(cè)器陣列為1040 mm，射線源到旋轉(zhuǎn)中心的距離為570 mm。重建的CT圖像大小為512像素×512像素，每個(gè)像素的大小為1.2 mm×1.2 mm。

設(shè)定不同的間隔掃描角度，從CT掃描中分別得獲得全投影(360個(gè)投影)、120個(gè)、90個(gè)和60個(gè)投影數(shù)的原始數(shù)據(jù)，利用FBP方法對(duì)這些數(shù)據(jù)重建得到對(duì)應(yīng)的重建圖像。全投影重建的圖像作為完整的圖像，用于和稀疏投影重建的圖像作對(duì)比。針對(duì)某一個(gè)投影數(shù)的重建圖像，從重建的圖像中隨機(jī)選取400張512*512的圖像作為訓(xùn)練數(shù)據(jù)，再從剩余中選取20張512*512的圖像作為測(cè)試數(shù)據(jù)。對(duì)于殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)的輸入，使用全投影圖像和稀疏重建后圖像之間的殘差作為訓(xùn)練初始化值，同時(shí)對(duì)不同投影數(shù)的重建圖像分別進(jìn)行訓(xùn)練，經(jīng)過20次訓(xùn)練后的殘差圖作為輸出。

2.2 網(wǎng)絡(luò)訓(xùn)練和參數(shù)的設(shè)定

本文為了訓(xùn)練的準(zhǔn)確性而得到較好的去偽影效果，不同投影數(shù)重建的訓(xùn)練樣本分別為400張，測(cè)試圖像分別為20張，同時(shí)把網(wǎng)絡(luò)的深度D設(shè)置為20。采用損失函數(shù)學(xué)習(xí)殘差映射μ(f)估計(jì)預(yù)測(cè)殘差n。優(yōu)化方法選擇Adam下降法，梯度的一階矩估計(jì)α為0.01，二階矩估計(jì)β1為0.9，β2為0.999，更新幅度的校正參數(shù)epsilon為1e-8。該殘差網(wǎng)絡(luò)的訓(xùn)練速率為10-3～10-4，每次訓(xùn)練逐漸遞減。為了產(chǎn)生較好的去偽影效果，提取出條狀偽影的特征圖，訓(xùn)練次數(shù)設(shè)置為60。

本文在訓(xùn)練的過程中使用MatConvNet工具包，運(yùn)行環(huán)境為MATLAB 2015a，處理器為Intel(R) Core(TM) i7_2600 CPU@ 3.4 GHz，內(nèi)存為8 GB，GPU顯卡為GeForce GTX 1070。在此配置環(huán)境中，通過該網(wǎng)絡(luò)訓(xùn)練樣本大約需要24個(gè)小時(shí)。對(duì)不同的投影數(shù)重建，通過測(cè)試程序恢復(fù)出重建的圖像，所需的測(cè)試時(shí)間見表1。

表1 不同投影數(shù)經(jīng)過稀疏重建恢復(fù)過程的測(cè)試時(shí)間

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

2.3.1 定性分析

分別利用120個(gè)投影、90個(gè)投影和60個(gè)投影重建的圖像作為樣本進(jìn)行訓(xùn)練，而后通過殘差圖像恢復(fù)清晰的CT圖。圖4顯示了不同投影重建圖像訓(xùn)練后恢復(fù)的結(jié)果，應(yīng)用本文提出的方法達(dá)到了很好的去偽影效果。

圖4 不同投影數(shù)重建圖像經(jīng)過殘差學(xué)習(xí)恢復(fù)的結(jié)果

隨著稀疏重建投影數(shù)的減少，仍能很好地恢復(fù)出重建的圖像，效果接近全投影的原圖像。使用殘差網(wǎng)絡(luò)訓(xùn)練方法達(dá)到了預(yù)期的處理目的了。

2.3.2 定量分析

為了驗(yàn)證本文提出的方法，觀察圖像像素強(qiáng)度尤為重要。本實(shí)驗(yàn)采用120個(gè)投影的重建圖像、90個(gè)投影的重建圖像和60個(gè)投影的重建圖像與全投影、恢復(fù)出的圖像作對(duì)比。選取全投影重建的圖像與對(duì)應(yīng)的不同投影個(gè)數(shù)重建的圖像和恢復(fù)出的圖像，測(cè)試其在同一幅CT圖像中200行，120列到250列像素的強(qiáng)度。

圖5中120個(gè)、90個(gè)和60個(gè)投影重建的圖像強(qiáng)度偏離了全投影的像素強(qiáng)度，尤其是60個(gè)投影重建的圖像偏離的比較嚴(yán)重，通過本文提出的方法恢復(fù)的圖像強(qiáng)度已經(jīng)很接近全投影的圖像強(qiáng)度了。

圖5 全投影、不同投影個(gè)數(shù)的重建和恢復(fù)的圖像(圖像矩陣的第200行)在不同像素位置的強(qiáng)度

2.3.3 圖像評(píng)價(jià)

在評(píng)價(jià)重建圖像的好壞主要有兩個(gè)指標(biāo)：峰值信噪比(peak signal to noise ratio，PSNR)和結(jié)構(gòu)相似性(structu-ral similarity，SSIM)[13]。

2.3.3.1 PSNR

PSNR的計(jì)算如下

(6)

(7)

MSE表示稀疏重建圖像X和通過殘差學(xué)習(xí)恢復(fù)的圖像Y的均方誤差(mean square error)，H、W分別為圖像的高度和寬度；n為像素的位深，在CT圖像中位深為12。

2.3.3.2 SSIM

SSIM是一種全參考的圖像質(zhì)量評(píng)價(jià)指標(biāo)，它分別從亮度、對(duì)比度、結(jié)構(gòu)3方面度量圖像相似性

SSIM(X,Y)=[l(X,Y)]α·[c(X,Y)]β·[s(X,Y)]γ

(8)

l(X,Y)為亮度對(duì)比函數(shù)，c(X,Y)為對(duì)比度對(duì)比函數(shù)，s(X,Y)為結(jié)構(gòu)對(duì)比度函數(shù)，X和Y分別表示稀疏重建的圖像和殘差學(xué)習(xí)恢復(fù)的圖像，α、β、γ是3個(gè)對(duì)比函數(shù)所占的比例因子。

本實(shí)驗(yàn)選取了20個(gè)原始數(shù)據(jù)作為測(cè)試數(shù)據(jù)，針對(duì)不同投影個(gè)數(shù)計(jì)算使用本文方法與全投影圖像比較的平均PSNR和SSIM值，結(jié)果見表2。

表2 不同投影數(shù)稀疏重建后使用殘差學(xué)習(xí)方法與全投影重建圖像比較的PSNR和SSIM平均值

表2反映的是訓(xùn)練恢復(fù)出的圖像與全投影重建后的圖像峰值信噪比和結(jié)構(gòu)相似度值。由上表的數(shù)據(jù)對(duì)比，隨著投影數(shù)的減少重建效果在一定程度上下降了。PSNR介于20 dB～40 dB之間，效果是最佳的。當(dāng)投影數(shù)太少，PSNR較小，效果要變得稍微差一些。SSIM的值也在逐漸下降。其中SSIM的值越接近于1，效果越好。

2.3.4 實(shí)驗(yàn)比較

為了驗(yàn)證本文提出方法的效果，實(shí)驗(yàn)部分采用了與傳統(tǒng)方法的比較。使用張俊峰等的伽瑪正則化[14]的CT稀疏重建方法分別對(duì)120個(gè)、90個(gè)和60個(gè)投影的數(shù)據(jù)進(jìn)行重建，觀察重建圖的效果和相應(yīng)的PSNR、SSIM值。圖6從左到右依次是稀疏投影重建的圖像、使用伽瑪正則化方法重建的圖像和本文提出的方法重建的圖像。

從測(cè)試數(shù)據(jù)的平均PSNR和SSIM分析，可以明顯的得出本文提出的方法要好于傳統(tǒng)的伽瑪正則化方法，見表3、表4。在一定程度上，利用殘差學(xué)習(xí)的方法去偽影較伽瑪正則化方法好。

表3 不同投影數(shù)的平均PSNR/dB比較

表4 不同投影數(shù)的平均SSIM比較

圖6 伽瑪正則化方法和本文方法的比較

3 結(jié)束語

與傳統(tǒng)方法不同的是本文研究的方法在稀疏重建的基礎(chǔ)上引入了殘差學(xué)習(xí)的結(jié)構(gòu)框架，在一定的程度上大大降低了X射線的輻射劑量。同時(shí)采用GPU加速，縮短了訓(xùn)練的時(shí)間，加速了實(shí)驗(yàn)的訓(xùn)練過程。

總之，本文提出的校正方法不需要大面積掃描人體，在現(xiàn)有的條件下不增加設(shè)備，不影響圖像的成像質(zhì)量；降低了輻射對(duì)人體不必要的傷害，達(dá)到了很好的效果。另外，該殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)比較單一，需要進(jìn)一步優(yōu)化和改進(jìn)，這將是本文下一步的重點(diǎn)工作。