基于加權字典學習方法的低劑量CT圖像重建

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

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

引言

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

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

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

1 加權字典學習重建算法

1.1 字典學習重建模型

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

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

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

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

基于字典學習方法的重建問題可以表示為帶有正則約束項的最優化問題：

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

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

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

1.2 權重因子

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

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

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

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

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

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

1.3 最優化問題求解步驟

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

迭代截止條件為：

2 實驗結果與分析

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

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

在未加噪聲的體模仿真實驗中，Shepp-Logan標準體模的掃描和重建參數，見表1，重建結果，見圖1；人體頭部切片的掃描和重建參數與Shepp-Logan實驗基本相同，不同之處是投影角度數變為180°，重建結果，見圖2。

表1 Shepp-Logan體模掃描與重建參數

圖1 Shepp-Logan標準體模仿真實驗重建結果（無噪聲）

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

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

表2 無噪聲數據重建圖像客觀評價

兩種字典學習方法在重建頭部切片時的收斂速率，見圖3。可以看到加權字典學習算法略微慢于未加權的算法，因為權重因子削弱了正則約束項對圖像細節區域的平滑效果，使其收斂速度減慢。但與提高的成像質量相比，微小的收斂速度下降是可以接受的。

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

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

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

表3 有噪聲數據重建圖像客觀評價

3 結論

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

[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] 焦鵬飛,李亮,趙驥.壓縮感知在醫學圖像重建中的最新進展[J].CT理論與應用研究,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.