基于正弦圖分區(qū)修復(fù)的稀疏角度CT重建算法

張 萌， 梁亞星， 陳 燕， 桂志國,， 張 權(quán),*

(1.中北大學電子測試技術(shù)國家重點實驗室， 太原 030051； 2.中北大學生物醫(yī)學成像與影像大數(shù)據(jù)山西省重點實驗室， 太原 030051； 3.中北大學信息與通信工程學院， 太原 030051)

計算機斷層(computed tomography, CT)成像技術(shù)在輔助臨床診斷方面具有巨大的應(yīng)用價值[1]。新冠狀病毒肺炎(COVID-19)防控期間，胸部CT檢查在患者的排查、診斷和預(yù)后評估等方面均發(fā)揮著不可替代的作用[2]。然而，常規(guī)CT的輻射劑量較高，其累積輻射會為患者帶來潛在的致癌風險[3]。因此，開展有關(guān)低劑量掃描技術(shù)的研究尤為重要。目前存在的諸多降低劑量的措施中，稀疏角度掃描是較重要的一種，通過減少投影采樣角度，在一定程度上降低了X射線的照射劑量。然而，該掃描方式造成的投影數(shù)據(jù)缺失，會導(dǎo)致重建圖像含有嚴重的條形偽影。

針對該問題，研究人員提出了多種稀疏重建算法，可概括為兩類：一類是迭代重建算法，如宋潔等[4]用LP范數(shù)正則項代替重建模型中的全變差(total variation, TV)正則項，實現(xiàn)了36個掃描角度下的精確CT重建。此類方法重建精度高，但計算量大且耗時。相比而言，基于投影數(shù)據(jù)恢復(fù)的算法先補全缺失的投影數(shù)據(jù)，再進行解析或迭代重建，重建速度較快。許多可行的研究皆采用此類方法實現(xiàn)快速重建[5]。Li等[6]對每個缺失投影數(shù)據(jù)均擬合出一組經(jīng)過該點的正弦曲線，并基于特征向量選取最佳曲線引導(dǎo)插值，提高了數(shù)據(jù)補全準確度并改善了重建圖像的結(jié)構(gòu)模糊問題；Kim等[7]使用全變差去噪算法對類正弦曲線進行分解，驗證了該方法用于投影正弦圖補全的可行性；Li等[8]采用字典學習方法修復(fù)稀疏投影數(shù)據(jù)，主要利用了投影數(shù)據(jù)的稀疏性。然而，該算法忽視了投影數(shù)據(jù)間的非局部相似性，易造成結(jié)構(gòu)細節(jié)模糊。近年來，深度學習技術(shù)的爆發(fā)式發(fā)展也使其在低劑量CT重建領(lǐng)域有了初步應(yīng)用[9]。Ghani等[10]及韋子權(quán)等[11]分別將深度學習網(wǎng)絡(luò)用于稀疏角度CT重建的投影域修復(fù)及圖像后處理。然而不容忽視的是醫(yī)療行業(yè)在倫理及信息安全的限制，使得獲取大規(guī)模醫(yī)療數(shù)據(jù)集較困難。相較而言，傳統(tǒng)算法的研究則限制較少，因此，對傳統(tǒng)稀疏重建算法的研究仍具價值。

鑒于此，在Li等[8]工作的基礎(chǔ)上，提出一種基于正弦圖分區(qū)修復(fù)的稀疏角度CT重建算法。分析正弦圖不同區(qū)域的特性；利用此特性對正弦圖進行分區(qū)稀疏表示修復(fù)；在稀疏表示模型中加入低秩懲罰，以便引入正弦圖的非局部信息，進而保留正弦圖的內(nèi)部結(jié)構(gòu)。

1 正弦圖特性分析

圖1(a)所示為常規(guī)掃描的臨床投影正弦圖，可以觀察到正弦圖邊界區(qū)域的投影值偏小，內(nèi)部區(qū)域投影值偏大且差異較小。若將正弦圖劃分為若干n×n的重疊子塊，則邊界子塊中多為零像素，故該區(qū)域圖像塊的稀疏性較強；在內(nèi)部區(qū)域隨機選取2 000個圖像塊，分別展開成列向量并組成新矩陣，繪出該矩陣的奇異值分布，如圖1(b)所示。易見，奇異值衰減極快，說明該矩陣具有低秩性，即正弦圖的內(nèi)部具有非局部相似性。

圖1 原始正弦圖及奇異值分布Fig.1 Original sinogram and the distribution of singular values

2 算法

為了充分利用不同區(qū)域的特性，對稀疏采樣的正弦圖進行分區(qū)修復(fù)。采用圖像灰度熵對正弦圖像塊分類并訓(xùn)練分類字典，利用字典學習算法修復(fù)邊界子塊，構(gòu)建新的聯(lián)合修復(fù)模型修復(fù)內(nèi)部子塊。內(nèi)部子塊是構(gòu)成正弦圖的主要部分，因而算法的關(guān)鍵在于該部分的修復(fù)。算法框圖如圖2所示。

圖2 算法框架Fig.2 The algorithm framework

2.1 基于灰度熵的圖像塊分類

采用灰度熵區(qū)分不同的正弦圖子塊，灰度熵較大的為邊界子塊；反之，則為內(nèi)部子塊。若大小為n×n圖像塊中像素的灰度級為n2，pj表示灰度級j出現(xiàn)的概率，則灰度熵h定義為

(1)

鑒于稀疏投影中有效數(shù)據(jù)少，增加線性插值操作，將插值后的圖像塊與稀疏采樣的正弦圖像塊對應(yīng)組成圖像塊對，利用前者的灰度熵對后者分類。

2.2 正弦圖邊界子塊修復(fù)

正弦圖邊界子塊稀疏性較強，對其采用字典學習方法進行修復(fù)[8]。

若存在字典D∈Rn2×p(R為實數(shù)域)，正弦圖像塊對應(yīng)向量u∈Rn2的稀疏表示為

u=Dα

(2)

式(2)中：α為稀疏表示向量。

若z∈Rm2為稀疏采樣的正弦圖像塊向量，則

z=Mu+ε

(3)

式(3)中：M∈Rm2×n2為二進制下采樣矩陣(標記了已知和缺失元素位置)；ε為噪聲。

在式(2)、式(3)的基礎(chǔ)上，基于字典學習的正弦圖像塊修復(fù)算法[8]可總結(jié)為

z=MDαz+ε

(4)

(5)

u=Dαz

(6)

其思想如下：原始正弦圖像塊與稀疏采樣正弦圖像塊具有近似的稀疏表示。使用易求的αz代替α與字典D相乘，即可修復(fù)正弦圖像塊。

2.3 正弦圖內(nèi)部子塊修復(fù)

為了充分利用正弦圖的非局部信息且提高稀疏表示性能，受文獻[12-13]啟發(fā)，在稀疏表示模型中引入了低秩矩陣填充模型，同時考慮了正弦圖的局部稀疏性與非局部低秩性。

在由一個不完整的矩陣恢復(fù)出一個完整的低秩矩陣過程中，標準矩陣填充模型的求解常被松弛為以下凸約束優(yōu)化模型[14]：

(7)

式(7)中：L表示已經(jīng)觀測到的不完整矩陣；Ω表示觀測值的索引集合；PΩ(·)為投影算子，表示將空間域Ω外的元素全部置0，Ω以內(nèi)的元素保持不變。

鑒于傳統(tǒng)低秩矩陣填充模型不能解決稀疏投影數(shù)據(jù)這樣的整行或整列元素缺失的問題[12]，將其融入稀疏表示模型中，通過構(gòu)建稀疏正弦圖低秩矩陣，整合正弦圖非局部信息，所得聯(lián)合修復(fù)模型為

s.t.X=DA,PΩ(X)=PΩ(Z)

(8)

式(8)中：Z=Rn2×r表示已知的稀疏采樣矩陣(缺失數(shù)據(jù)為0)；A∈Rp×r為稀疏系數(shù)矩陣；β為懲罰參數(shù)。

構(gòu)建低秩矩陣Z時，無需搜索相似塊，直接將稀疏采樣正弦圖內(nèi)部子塊各自展開成列向量，再組成新矩陣。但這樣得到的低秩矩陣列數(shù)過于龐大，因此，采取分批處理的方法。

2.4 模型求解

對于上述模型，采用交替方向法求解，將式(8)轉(zhuǎn)化為交替求解兩個子問題：

(9)

(10)

式(9)為一個稀疏表示問題，利用正交匹配追蹤(OMP)算法求解[15]。式(10)為一個矩陣填充問題，采用增廣拉格朗日乘子法(ALM)求解[12]，其增廣拉格朗日函數(shù)表示為

(11)

式(11)中：μ為拉格朗日乘子。

最小化上述增廣拉格朗日函數(shù)，可以轉(zhuǎn)換成交替更新如下問題：

(12)

式(12)中：X的更新可以表示為

(13)

式(13)可以使用奇異值閾值(SVT)算法進行求解，得到

(14)

式(14)中：Wk=DA-Yk/μk；Uk、Vk分別為Wk的左、右奇異矩陣，即Wk=UkSk(Vk)T。

求解X問題時，Ω空間域內(nèi)元素為已知測量值，即PΩ(X)=PΩ(Z)，因此X的最終更新方式為

(15)

2.5 算法描述

算法具體步驟概括如下。

(1)計算訓(xùn)練集中圖像塊的灰度熵，依據(jù)灰度熵大小將訓(xùn)練集分為兩類，分別訓(xùn)練邊界字典Db與內(nèi)部字典Ds。

(2)進行預(yù)處理并分類。對稀疏采樣數(shù)據(jù)進行線性插值，將插值前后的圖像分別劃分為n×n的重疊子塊并剔除無效塊(像素值全為0)，根據(jù)插值后圖像塊的灰度熵值，對待修復(fù)的正弦圖子塊進行分類。

(3)正弦圖邊界子塊修復(fù)。

輸入：字典Db，下采樣矩陣M，含缺失元素的邊界子塊。

針對每一子塊，依次計算式(5)和式(6)。

輸出：已修復(fù)的正弦圖子塊。

(4)正弦圖內(nèi)部子塊修復(fù)。利用待修復(fù)的內(nèi)部正弦圖子塊構(gòu)建低秩矩陣Z，并將其劃分為T個批次，即Z={Z1,Z2,…,ZT}，針對每一子矩陣Zi，采用算法1稀疏-低秩聯(lián)合修復(fù)算法補全缺失像素值。

(5)將步驟(3)、(4)中已修復(fù)的圖像子塊及低秩矩陣，映射到正弦圖中對應(yīng)的位置，重疊區(qū)域取平均處理，即可得到補全的完整正弦圖，最后進行FBP重建，得到重建圖像。

算法1稀疏-低秩聯(lián)合修復(fù)算法

輸入：字典Ds，下采樣矩陣M，構(gòu)建的低秩子矩陣Zi。

輸出：已修復(fù)的低秩矩陣X。

初始化：X0=Zi，計算式(5)將結(jié)果作為A的初值，Y0=0，μ0>0，ρ>1，β>0，k=1

whilek

求解式(10)，得到Xk+1

求解式(9)，得到Ak+1

k=k+1

end while

3 實驗結(jié)果

分別采用模擬頭模數(shù)據(jù)與實際骨盆數(shù)據(jù)進行實驗驗證。算法運行壞境如下：CPU為Intel(R) Core(TM)i5-4210H；內(nèi)存為12 G；操作系統(tǒng)為Windows 10；編程平臺為Matlab R2016a。

采用線性插值算法、字典學習算法[8]，基于圖像灰度熵的自適應(yīng)字典學習算法[16]作為對比算法。使用完整的投影數(shù)據(jù)構(gòu)建訓(xùn)練集，并采用K-奇異值分解(K-SVD)算法[17]訓(xùn)練字典。修復(fù)的投影數(shù)據(jù)采用基于ramp濾波器的FBP算法重建。

實驗中，采用均方誤差(MSE)、峰值信噪比(PSNR)與視覺信息保真度(VIF)3種量化指標對修復(fù)的正弦圖及重建圖像進行定量分析。

3.1 模擬數(shù)據(jù)實驗

選擇Shepp-Logan頭模模擬生成平行束投影數(shù)據(jù)，頭模大小為256像素×256像素，探測單元的數(shù)量為260個。對完整投影數(shù)據(jù)每間隔4°采集一個樣本(180°內(nèi)均勻采集45個角度)模擬稀疏掃描。實驗中的相關(guān)參數(shù)設(shè)置如下。

字典參數(shù)：圖像塊大小=8×8；字典大小=64×256；邊界稀疏度=3；內(nèi)部稀疏度=4；迭代次數(shù)為10。

其他參數(shù)：β=1；c=10；ρ=1.1；T=10；灰度熵閾值t=0；迭代次數(shù)N=50；μ0=c/max[δ(Z)]，其中δ(Z)表示矩陣Z的奇異值。

圖3(a)為稀疏采樣正弦圖，圖3(b)～圖3(f)為原始完整正弦圖和4種算法修復(fù)的正弦圖及相應(yīng)的感興趣區(qū)域放大圖。觀察放大圖可以發(fā)現(xiàn)：圖3(c)邊緣結(jié)構(gòu)處呈現(xiàn)出鋸齒狀；圖3(d)存在明顯的豎直條紋；與圖3(d)相比圖3(e)的視覺效果有所提高，但在部分區(qū)域仍含有豎直條紋；相比而言，本文算法修復(fù)的正弦圖在邊緣與平坦區(qū)域均呈現(xiàn)較好的一致性。

此外，在圖4中繪制了各修復(fù)正弦圖在第140行的剖面輪廓，可以看出所提算法修補的正弦圖與原始正弦圖的輪廓最為吻合，表明本文算法具有較好的正弦圖修復(fù)效果。

圖5(a)為圖3(a)的FBP重建結(jié)果，易見，稀疏重建圖像中引入了大量的條形偽影，圖像退化嚴重。圖5(b)～圖5(f)為圖3中原始正弦圖及各修復(fù)正弦圖的FBP重建圖像及其局部放大圖。可見，圖5(c)中依然包含清晰的條形偽影，且結(jié)構(gòu)模糊較嚴重；與前者相比，圖5(d)與圖5(e)中的條形偽影及結(jié)構(gòu)模糊問題均有改善，但箭頭標記處均含有環(huán)形偽影；相比而言，本文算法的重建結(jié)果偽影較少且矩形標注區(qū)域結(jié)構(gòu)的清晰度有所提高。

表1列出了各算法修復(fù)正弦圖所用時間及質(zhì)量參數(shù)對比。相應(yīng)的在表2中列出了各修復(fù)正弦圖經(jīng)FBP重建后的圖像質(zhì)量參數(shù)對比。綜合來看，4種算法中本文算法雖耗時最長，但獲得的PSNR與VIF最高，MSE最小，算法的有效性得到驗證。

圖3 采用不同算法修復(fù)的正弦圖及其感興趣區(qū)域放大圖Fig.3 Repaired sinograms and their regions of interest enlargement images by different methods

圖4 各修復(fù)正弦圖第140行的剖面輪廓對比Fig.4 Profile contour contrast of the repaired sinograms along the 140th row

表1 模擬數(shù)據(jù)實驗中正弦圖修復(fù)時間及質(zhì)量參數(shù)對比Table 1 Comparison of repair time and quality parameters of sinogram in simulated data experiment

圖5 圖3中已修復(fù)正弦圖的重建結(jié)果及局部放大圖Fig.5 Reconstruction results and their partial enlargement images of the repaired sinograms in Fig.3

表2 模擬數(shù)據(jù)實驗中重建圖像的質(zhì)量參數(shù)對比Table 2 Quality parameters comparison of reconstructed images in simulated data experiment

3.2 實際數(shù)據(jù)實驗

采用扇束掃描實際骨盆投影數(shù)據(jù)為實驗對象，重建圖像大小為256像素×256像素, 探測單元的數(shù)量為512個。對完整正弦數(shù)據(jù)每間隔4°采集一個樣本(360o內(nèi)均勻采集90個角度)，從而獲得稀疏采樣正弦圖。實驗的相關(guān)參數(shù)設(shè)置如下。

字典參數(shù)：圖像塊大小為11×11；字典大小為121×484；其余同模擬數(shù)據(jù)。

其他參數(shù)：T=36；迭代次數(shù)N=60；其余同模擬數(shù)據(jù)。

圖6(a)為實際數(shù)據(jù)的稀疏采樣正弦圖，圖6(b)～圖6(f)為原始完整正弦圖和4種算法修復(fù)的正弦圖及其感興趣區(qū)域放大圖。對比矩形標注區(qū)域，本文算法在邊緣結(jié)構(gòu)與圖像一致性的保持兩方面均有較好表現(xiàn)，其修復(fù)的正弦圖[圖6(f)]的視覺效果優(yōu)于另3種算法。

圖7為各修復(fù)正弦圖在第340行的剖面輪廓對比圖，為方便觀察，僅顯示第250列到第360列。顯然，本文算法修復(fù)的正弦圖不論在邊緣跳變部分還是平坦區(qū)域皆與原始正弦圖的接近程度最高。

圖6 采用不同算法修復(fù)的正弦圖及其局部放大圖Fig.6 Repaired sinograms and their partial enlargement images by different algorithms

圖7 各修復(fù)正弦圖第340行的剖面輪廓對比Fig.7 Profile contour contrast of the repaired sinograms along the 340th row

圖8(a)為稀疏投影的FBP重建結(jié)果，可見圖中含有嚴重的條形偽影。圖8(b)～圖8(f)給出了圖6中原始正弦圖與各修復(fù)正弦圖的FBP重建結(jié)果，同時選取圖中矩形標注區(qū)域放大顯示。可以發(fā)現(xiàn)，線性插值算法修復(fù)的正弦圖中條形偽影去除不足，兩種字典學習算法修復(fù)的正弦圖皆存在過平滑現(xiàn)象，正弦圖的部分結(jié)構(gòu)細節(jié)丟失。較之于前3種算法，本文算法有效抑制了條形偽影且改善了結(jié)構(gòu)模糊與過平滑問題，在稀疏角度重建中更具優(yōu)勢。更進一步，對4種重建算法進行定量分析，如表3和表4所示，對比修復(fù)時間及PSNR、MSE和VIF的變化亦可得到相同結(jié)論，算法的有效性得到進一步驗證。

圖8 圖6中已修復(fù)正弦圖的重建結(jié)果及局部放大圖Fig.8 Reconstruction results and their partial enlargement images of the repaired sinograms in Fig.6

表3 實際數(shù)據(jù)實驗中正弦圖修復(fù)時間及質(zhì)量參數(shù)對比Table 3 Comparison of repair time and quality parameters of sinogram in real data experiment

表4 實際數(shù)據(jù)實驗中重建圖像的質(zhì)量參數(shù)對比Table 4 Quality parameters comparison of reconstructed images in real data experiment

4 結(jié)論

區(qū)別于傳統(tǒng)圖像稀疏表示只訓(xùn)練單一字典，本文算法依據(jù)各正弦圖子塊的不同特性，分別訓(xùn)練兩個字典，提高了各類子塊的稀疏表示能力。修復(fù)內(nèi)部子塊時，構(gòu)建了聯(lián)合修復(fù)模型，該模型同時考慮了正弦數(shù)據(jù)間的局部稀疏性與非局部低秩性，可以進一步保留正弦圖的內(nèi)部結(jié)構(gòu)。仿真及臨床數(shù)據(jù)實驗表明，較之于3種傳統(tǒng)的重建算法，本文算法具有以下特點：

(1)正弦圖的修復(fù)能力有所提高，其邊緣及平坦區(qū)域的一致性更強。

(2)抑制稀疏角度CT重建圖像條形偽影的同時改善了結(jié)構(gòu)模糊及過度平滑問題。

然而，與其他方法相比，本文算法在修復(fù)正弦圖時需要更長的時間，下一步工作將考慮采用并行化方法加速修復(fù)。