基于加權(quán)字典學(xué)習(xí)方法的低劑量CT圖像重建

章程，張健，杜強，李銘，劉景鑫

1.中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇 蘇州 215163；2.長春市計量檢定測試技術(shù)研究院，吉林 長春 130012；3.吉林大學(xué)中日聯(lián)誼醫(yī)院 放射線科，吉林 長春 130033

引言

由于計算機斷層成像（Computed Tomography，CT）空間分辨率高、圖像質(zhì)量好，是一種廣泛使用的醫(yī)療檢查方法。然而，X射線的過度輻射會導(dǎo)致基因疾病和癌癥[1]。降低輻射劑量主要通過降低X射線輻射強度或者減少采樣率來實現(xiàn)，前者會造成重建圖像噪聲增加，信噪比降低，而后者是一個欠采樣不完備投影數(shù)據(jù)重建問題。當(dāng)利用以濾波反投影[2]算法為代表的解析重建算法進行重建時，由于采樣率低于香濃奈奎斯特采樣定理[3]的要求，重建圖像中會有嚴(yán)重的混疊偽影，極大的影響了成像質(zhì)量。

應(yīng)對這種欠采樣的重建問題，迭代重建算法有著極大的優(yōu)勢。2006年，由Candes等[4-5]提出的壓縮感知（Compressed Sensing，CS）理論開拓了迭代算法的新方向[6]。該理論表明，信號可以在特定變換域被稀疏表示，從而可以突破采樣定理的限制，用較低的采樣頻率觀測和重構(gòu)信號。以CS理論為基礎(chǔ)，如今已出現(xiàn)了一些基于稀疏約束的欠采樣重建算法。其中一類以全變分（Total Variation，TV）最小化為稀疏約束，具有代表性的算法有自適應(yīng)最速下降—凸集投影（Adaptive Steep est Descent-Projection onto Convex Sets，ASD-POCS）算法[7]和梯度投影Barzilai-Borwein（Gradient Projection Barzilai-Borwein，GPBB）[8]算法。另一類算法的稀疏正則約束項由字典學(xué)習(xí)（Dictionary Learning，DL）方法生成，該方法是將待重建的CT圖像分解為互相重疊的小圖像塊，這些圖像塊都可以用一個過完備字典（列數(shù)遠(yuǎn)大于行數(shù)的矩陣，每一列都叫做這個字典的原子）稀疏表示。2012年，Xu等[9]和Bai等[10]提出了基于字典學(xué)習(xí)方法的低劑量CT重建算法，得到了不錯的重建結(jié)果。

然而，這些典型的稀疏重建算法在利用正則約束項抑制噪聲的同時，在圖像的低對比度信息區(qū)域也引入了相同的平滑效果。這種過平滑效應(yīng)使圖像的空間分辨率降低，細(xì)節(jié)信息丟失。為了保留更多的邊緣細(xì)節(jié)信息，基于TV最小化的算法通過加入合理的自適應(yīng)權(quán)重因子使算法能夠自動分辨噪聲和圖像細(xì)節(jié)信息，從而抑制圖像中低對比度信息區(qū)域的平滑效果，保留軟組織邊緣信息[11-12]。本文采用類似的優(yōu)化策略，在字典學(xué)習(xí)算法模型中加入權(quán)重因子。通過對比Shepp-Logan投影數(shù)據(jù)和人體頭部切片的重建實驗結(jié)果，與未采用權(quán)重因子的算法相比，本文提出的加權(quán)字典學(xué)習(xí)算法重建質(zhì)量更好，空間分辨率更高，噪聲抑制效果更好。與聯(lián)合代數(shù)重建（Simultaneous Algebraic Reconstruction Technique，SART）算法[13]和GPBB算法的結(jié)果對比也體現(xiàn)出了加權(quán)字典學(xué)習(xí)算法的優(yōu)越性。

1 加權(quán)字典學(xué)習(xí)重建算法

1.1 字典學(xué)習(xí)重建模型

CT圖像 的重建問題是一個典型的逆問題，投影數(shù)據(jù)采集過程可以表示為一個線性方程

其中 A={aij}∈RI×N2為系統(tǒng)矩陣；μ∈RN2×1為待求解的線性衰減系數(shù)分布，也就是待重建的CT圖像，圖像大小為N×N像素，重排得到一維向量μ；∈RI×1為測量的線積分投影數(shù)據(jù)。

在單色X射線源的假設(shè)下，CT投影測量值符合泊松分布規(guī)律：

其中b∈RI×1和 y∈RI×1分別為入射前和透射后的 X 射線強度；g=(g1,g2,K,gI)T∈RI×1為真實的線積分投影數(shù)據(jù)，其測量值為= ln(bi/(yi? γi))；γi為讀出噪聲。

基于字典學(xué)習(xí)方法的重建問題可以表示為帶有正則約束項的最優(yōu)化問題：

其中ωi=(yi?γi)2/yi為統(tǒng)計迭代模型中的統(tǒng)計權(quán)重；是從N×N維圖像中抽取N0×N0維小圖像塊的算子；圖像塊相互重疊，當(dāng)抽取算子的滑動距離為1時，圖像塊的總數(shù)S=(N-N0+1)×(N-N0+1)；為過完備字典被稱為字典中的原子；as∈RK×1是圖像塊以字典D為基底的稀疏表示；λ與vs為正則化參數(shù)。最優(yōu)化問題利用交替最小化方法求解，每一次迭代都分為兩個連續(xù)的步驟——字典更新與圖像更新。字典更新過程中，圖像保持不變，目標(biāo)函數(shù)簡化為：

分別利用K-SVD（K Singular Value Decomposition）[14]算法和正交匹配追蹤（Orhtogonal Matching Pursuit，OMP）[15]算法更新字典D和稀疏表示as，OMP算法設(shè)定稀疏度為一個定值，比如5～10，將式(4)的L0范數(shù)約束優(yōu)化問題簡化為OMP問題，此時求解不需要設(shè)定正則化參數(shù)vs的值。圖像更新過程中，字典和稀疏表示不變，目標(biāo)函數(shù)簡化為：

利用可分離二次代理函數(shù)[16]方法，圖像更新公式可以表示為：

1.2 權(quán)重因子

在字典學(xué)習(xí)重建模型中，式(3)中的正則約束項為：

上式中，對于不同的圖像塊，正則約束項的權(quán)重都相同，因此圖像在迭代更新的過程中算法對不同區(qū)域的平滑效果相同，無法有效的在抑制噪聲的同時保留圖像細(xì)節(jié)信息。為了提高重建圖像的質(zhì)量和空間分辨率，本文提出根據(jù)圖像塊中所含邊緣細(xì)節(jié)信息的多少選擇適當(dāng)?shù)臋?quán)重因子給予不同程度的平滑效果，圖像塊所含細(xì)節(jié)信息越多，權(quán)重因子越小。

由于圖像塊稀疏表示的求解通過OMP算法實現(xiàn)，且稀疏度為一個定值，因此圖像塊與其稀疏表示的殘差可以間接反映細(xì)節(jié)信息的多少，殘差越大，細(xì)節(jié)信息越多。本文選擇上一次迭代結(jié)果的殘差的反比例函數(shù)作為這一次迭代的權(quán)重因子，計算表達(dá)式為：

其中ε為保證算法穩(wěn)定性的小量，Ct是所有圖像塊稀疏表示殘差的平均值，是一個平衡參數(shù)，其作用是使當(dāng)前計算出的權(quán)重因子不影響到目標(biāo)函數(shù)中正則化參數(shù)λ的取值，加權(quán)字典學(xué)習(xí)重建模型的正則約束項變化為：

在每次迭代之前，權(quán)重因子根據(jù)上一次迭代結(jié)果自動更新，細(xì)節(jié)信息越多的圖像塊會傾向于得到一個越小的權(quán)重，使圖像的低對比度信息得到更好的保留。由于權(quán)重因子的加入，圖像更新過程的目標(biāo)函數(shù)，即式(5)變化為：

圖像更新公式(6)變化為：

1.3 最優(yōu)化問題求解步驟

本文提出的算法以圖像塊殘差作為權(quán)重因子的計算依據(jù)，稱為加權(quán)字典學(xué)習(xí)（Weighted Dictionary Learning，WDL）算法。本文算法用有序子集凸優(yōu)化[17]（Ordered Subsets Convex，OSC）算法進行迭代加速，其具體求解步驟如下所示：① 步驟1：初始化，μ0為隨機矩陣，D0為離散余弦變換矩陣，迭代次數(shù)t=0，權(quán)重因子=1；② 步驟2：用OSC算法迭代加速；③ 步驟3：用K-SVD算法更新字典得到Dt+1；④ 步驟4：用OMP算法更新稀疏表示得到；⑤ 步驟5：利用式(8)更新權(quán)重因子得到；⑥ 步驟6：利用式(11)更新圖像，得到μt+1，迭代次數(shù)t=t+1；⑦ 步驟7：若滿足迭代截止條件，停止迭代，返回重建結(jié)果，否則，返回步驟2。

迭代截止條件為：

2 實驗結(jié)果與分析

為了評估算法的性能，本文分別用Shepp-Logan標(biāo)準(zhǔn)體模和真實人體頭部切片做數(shù)據(jù)仿真實驗，投影幾何為扇形束投影，投影角度均勻分布在360°范圍內(nèi)。對比算法為未加權(quán)重的字典學(xué)習(xí)算法，GPBB算法和SART算法。圖像質(zhì)量的客觀評價標(biāo)準(zhǔn)為歸一化平均絕對偏差（Normalized Mean Absolute Deviation，NMAD） 以 及 信 噪 比SNR：

可以看出，絕對偏差越小，信噪比越高，重建質(zhì)量越好。正則化參數(shù)λ的選取利用文獻[18]的自動選取模型。

在未加噪聲的體模仿真實驗中，Shepp-Logan標(biāo)準(zhǔn)體模的掃描和重建參數(shù)，見表1，重建結(jié)果，見圖1；人體頭部切片的掃描和重建參數(shù)與Shepp-Logan實驗基本相同，不同之處是投影角度數(shù)變?yōu)?80°，重建結(jié)果，見圖2。

表1 Shepp-Logan體模掃描與重建參數(shù)

圖1 Shepp-Logan標(biāo)準(zhǔn)體模仿真實驗重建結(jié)果（無噪聲）

通過對比發(fā)現(xiàn)，SART算法僅通過反投影差值平方最小化來重建圖像，沒有其他的先驗知識和約束條件，因此重建結(jié)果最差，含有很多偽影。GPBB算法是一種TV最小化算法，能夠出色的重建出圖像的平滑區(qū)域，但從差值圖像中可以看出在重建圖像的邊緣區(qū)域有過平滑效應(yīng)，因此差值圖像中有比較明顯的邊緣痕跡。未加權(quán)重的字典學(xué)習(xí)算法與GPBB算法相比，在重建上沒有優(yōu)勢，平滑效果不如GPBB算法，在圖像平滑區(qū)域有部分散粒噪聲，但字典學(xué)習(xí)算法在圖像邊緣區(qū)域的過平滑效應(yīng)也不如GPBB算法那么嚴(yán)重。本文提出的加權(quán)字典學(xué)習(xí)重建算法因為權(quán)重因子的引入，能夠自動區(qū)分重建圖像的平滑區(qū)域和邊緣細(xì)節(jié)區(qū)域，從而給予不同程度的平滑效果，其與真實圖像的差值圖像，既沒有明顯的散粒噪聲，也沒有邊緣痕跡，提高了重建圖像的空間分辨率和成像質(zhì)量。相比于未加權(quán)的字典學(xué)習(xí)算法，本文算法得到的結(jié)果的客觀評價標(biāo)準(zhǔn)也有了很大的提高，其結(jié)果，見表2。

圖2 真實人體頭部切片仿真實驗重建結(jié)果（無噪聲）

表2 無噪聲數(shù)據(jù)重建圖像客觀評價

兩種字典學(xué)習(xí)方法在重建頭部切片時的收斂速率，見圖3?？梢钥吹郊訖?quán)字典學(xué)習(xí)算法略微慢于未加權(quán)的算法，因為權(quán)重因子削弱了正則約束項對圖像細(xì)節(jié)區(qū)域的平滑效果，使其收斂速度減慢。但與提高的成像質(zhì)量相比，微小的收斂速度下降是可以接受的。

圖3 頭部切片重建NMAD隨迭代次數(shù)下降曲線圖

為了檢驗算法對噪聲的魯棒性，在人體頭部切片的仿真實驗中加入泊松噪聲，設(shè)定入射前的X射線強度為120萬光子數(shù)，透射物體后的探測光子數(shù)利用泊松分布模擬泊松噪聲。重建結(jié)果，見圖4?？陀^評價標(biāo)準(zhǔn)結(jié)果，見表3。結(jié)果與無噪聲的方針實驗結(jié)果相似，WDL算法結(jié)果的NMAD和SNR指標(biāo)與GPBB算法接近。但是通過觀察圖中黃色邊框內(nèi)的放大區(qū)域，發(fā)現(xiàn)GPBB算法的過平滑效應(yīng)導(dǎo)致圖像中有比較明顯的階梯效應(yīng)，其空間分辨率有所降低。從圖4g可以發(fā)現(xiàn)，DL算法在有噪聲的情況下對圖像的高對比度邊緣有過強的平滑效果。WDL算法雖然在一定程度上減弱了這種效果，但并沒有完美的解決這個問題，這也是WDL算法結(jié)果的評價標(biāo)準(zhǔn)略差于GPBB算法的原因。

圖4 真實人體頭部切片仿真實驗重建結(jié)果（有噪聲）

表3 有噪聲數(shù)據(jù)重建圖像客觀評價

3 結(jié)論

本文針對低劑量CT重建問題中，字典學(xué)習(xí)重建算法無法有效區(qū)分噪聲和低對比度信息，重建圖像容易丟失軟組織邊緣細(xì)節(jié)信息的問題，提出了一種通過權(quán)重因子保留細(xì)節(jié)信息，提高成像質(zhì)量的加權(quán)字典學(xué)習(xí)算法。權(quán)重因子與圖像塊用字典稀疏表示的殘差成反比，能夠有效抑制迭代過程中對圖像邊緣的過平滑效應(yīng)，同時保留了抑制噪聲的能力，使抑制噪聲和保留圖像細(xì)節(jié)信息不再是相互制約的關(guān)系。仿真實驗表明，提出的加權(quán)字典學(xué)習(xí)算法提升了字典學(xué)習(xí)算法的重建效果，其收斂速率相比未加權(quán)的算法僅有微小的下降。同時，與其他兩種對比算法相比，本文提出的算法在提高圖像空間分辨率上也有較大的優(yōu)勢。下一步研究需要探討DL算法在有噪聲的情形下對圖像高對比度邊緣的抑制效應(yīng)并對其進行改進。

[1] De gonzález AB,Mahesh M,Kim K,et al.Projected cancer risks from computed tomographic scans performed in the United States in 2007[J].Arch Intern Med,2009,169(22):2071-2077.

[2] Liao HY,Sapiro G.Sparse representations for limited data tomography[A].IEEE International Symposium on Biomedical Imaging: From Nano To Ma cro, IEEE[C].2008:1375-1378.

[3] Jerri AJ.The Shannon sampling theorem—Its various extensions and applications: a tutorial review[J].P IEEE,1977,65(11):1565-1596.

[4] Candes EJ,Romberg J,Tao T.Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency information[M].New York:IEEE Press,2006.

[5] Ca ndes EJ,Tao T.Near-optimal signal recovery from random projections: universal encoding strategies?[J].IEE E T Inform Theory,2006,52(12):5406-5425.

[6] 焦鵬飛,李亮,趙驥.壓縮感知在醫(yī)學(xué)圖像重建中的最新進展[J].CT理論與應(yīng)用研究,2012,21(1):133-147.

[7] Sidky EY,Pan XC.Image reconstruction in circular conebeam computed tomography by c onstrained, total-variation minimization[J].Phys Med Biol,2008,53(17):4777-4807.

[8] Justin CP,Song B,Kim JS ,et al.Fast compressed sensing-based CBCT r econstruction using Barzilai-Borwein formulation for application to on-line IGRT[J].Med Phys,201 2,39(3):1207-1217.

[9] Xu Q,Yu H,Mou X,et al.Low-dose X-ray CT reconstruction via dictionary learning[J].IEEE T Med Imaging,2012,31(9):1682-1697.

[10] Bai T,Yan H,Shi F,et al.WE-G-18A-04: 3D dictionary learning based statistical iterative reconstruction for low-dose cone beam CT imaging[J].Med Phys,2014,41(6):527.

[11] Tian Z,Jia X,Yu an K,et al.Low-dose CT reconstr uction via edge-preserving total variation regularization[J].Phys Med Biol,2011,56(18) :5949-5967.

[12] Guo W,Yin W.EdgeCS: edge guided compressive sensing reconstruction[A].Visual Communications and Image Processing. International Society for Optics and Photonics[C].201 0:77440L-77440L-10.

[13] Andersen AH,Kak AC.Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algo rithm[J].Ultrasonic Imaging,1984,6(1):81-94.

[14] Aharon M,Elad M,Bruckstein A.rmK-SVD: an algorithm for designing overcomplete dictionaries for sparse representation[J].IEEE T Signal Proces,2006,54(11):4311-4322.

[15] Tropp JA,Gil bert AC.Signal recovery from random measurements via orthogonal matching pursuit[J].IEEE T Inform Theory,2007,53(12):4655-4666.

[1 6] Elbakri IA,Fessler JA.Statistical image reconstruction for polyenergetic X-ray computed tomography[J].IEEE T Med Imaging,2002,21(2):89-99.

[17] Kamphuis C,Bee kman FJ.Accelerated iterative transmission C T reconstruction using an ordered subsets convex algorithm[J].IEEE T Med Imaging,1998,17(6):1101-1105 .

[18] Zhang C,Zhang T,Zheng J,et al.A Model of regularization parameter determination in low-dose X-Ray CT reconstruction based on dictionary learning[J].Comput Math Methods Med,2015:1-12.