張瑩芳，陳 平，潘晉孝

(1.中北大學(xué) 信息與通信工程學(xué)院，太原 030051；2.中北大學(xué) 信息探測與處理山西省重點實驗室，太原 030051)

0 引言

X射線計算機斷層成像(Computed Tomography，CT)是用X射線對物體進行斷層掃描，利用從多個角度獲得的投影數(shù)據(jù)重建物體截面信息的成像技術(shù)[1]。該技術(shù)以其快速采集和高空間分辨率等優(yōu)點在臨床診斷中得到了廣泛的應(yīng)用[2]。但CT掃描會產(chǎn)生大量的X射線輻射，過多的X射線輻射可能導(dǎo)致癌癥和遺傳性疾病[3]，為了降低CT掃描導(dǎo)致的患病風(fēng)險，稀疏角度CT成像技術(shù)越來越受到人們的關(guān)注[4]。

稀疏角度CT重建是一種通過減少投影角度的數(shù)量來降低輻射劑量的方法[5]，但受經(jīng)典的奈奎斯特采樣定理要求的限制，認為從欠采樣的數(shù)據(jù)重建斷層圖像是不可行的，現(xiàn)在得益于壓縮感知(Compressed Sensing, CS)理論的存在，許多不適定的逆問題包括稀疏角度CT問題都可以通過壓縮感知技術(shù)得到有效地解決[6]。但商業(yè)中應(yīng)用最廣泛的濾波反投影(Filtered Back-Projection, FBP)解析重建算法對投影數(shù)據(jù)的完備性有特定的要求，需要較高的采樣率來獲得滿意的圖像質(zhì)量[7]。相比而言，迭代重建算法對數(shù)據(jù)的完備性要求不高，可以通過引入圖像先驗約束信息，提高重建質(zhì)量[8]。

2006年，Sidky等引入全變分(TV)最小化的概念，提出了一種TVM-POCS(Total Variation Minimization-Projection Onto Convex Sets)算法，用于從稀疏采樣或稀疏角度投影數(shù)據(jù)進行CT圖像重建[9]。但由于圖像的分段常數(shù)假設(shè)，在圖像重建過程中，會產(chǎn)生階梯偽影以及圖像邊緣過度平滑的結(jié)果。因此，許多TV算法的變體被提出且被用于CT圖像重建。Kim H等提出了一種新的非局部全變分算子(NLTV)[10]；Jin等提出了一種基于各向異性全變分(ATV)的有限角CT圖像重建方法[11]；Wang等提出了一種新的迭代加權(quán)各向異性全變分方法(RWATV)[12]；Guo等提出了一個加權(quán)方向全變分(WDTV)最小化模型[13]；Qu等提出了一種新的基于梯度方向全變分(BDTV)最小化算法，自適應(yīng)選擇變化梯度方向計算方向差算子[14];Xi等提出了一種高階全變分(HOTV)重構(gòu)模型以提高其重建去噪效果，且在一定程度上保護圖像的邊緣信息[15];Kim等提出了一類新的非局部全變分，結(jié)合一階和二階導(dǎo)數(shù)，以很好地保持圖像的平滑度，平衡一階導(dǎo)數(shù)引起的階梯偽影現(xiàn)象[16]；Sagheer等提出了一種通過低秩張量建模和全變分正則化的方法對低劑量CT圖像進行去噪[17]；Li等提出了一種結(jié)合使用廣義全變分(PWLS-TGV)和字典學(xué)習(xí)(DL)的懲罰加權(quán)最小二乘法來恢復(fù)圖像[18]；Kai等將一種基于圖像貼片匹配的非局部三維剪切稀疏正則化用于稀疏角度CT圖像重建，有效抑制噪聲，提高重建圖像質(zhì)量[19]；時鴻雁等提出了一種基于導(dǎo)引圖像濾波(GIF)和截斷全變分的CT重建算法[20]。這些方法與傳統(tǒng)的TV模型相比，在一定程度上，有效抑制了噪聲，但在邊緣保護、圖像局部細節(jié)保護方面仍有所欠缺。

l0范數(shù)的高稀疏性，使得l0范數(shù)梯度正則化在圖像去噪和保護邊緣特征方面有著明顯優(yōu)勢。Xu等提出了一個有限角度X射線CT的重建模型，交替保邊延拓和保邊光滑算法(Alternating Edge-preserving Diffusion and Smoothing, AEDS)，包含兩個正則化項，在兩個方向起到保邊擴散和平滑的作用，并提出了一種交替最小化算法對模型進行近似求解，該方法在有限角度的重建過程中能較好地去除偽影，保留圖像邊緣，得到更好的重建圖像[21]。

針對在稀疏角度CT重建中僅僅使用一個約束，不能很好地兼顧去噪和保邊的問題，本文基于迭代重建算法，將AEDS應(yīng)用于稀疏角度的重建算法，對模型添加兩個正則項約束，提出了一種基于TV和AEDS的雙約束稀疏角度CT重建算法。該重建算法首先使用ART迭代重建算法，得到初步的迭代重建圖像，然后使用TV繼續(xù)重建圖像，最后使用AEDS，分別在x和y方向引入一個正則化項，進行一定的稀疏變換，達到邊緣保護的作用。

1 CT重建算法

1.1 CT重建模型

在理想的CT成像系統(tǒng)中，圖像重建滿足以下方程：

Au=p

(1)

式中，u為待重建圖像矢量；p為探測器在不同角度的投影數(shù)據(jù)；投影矩陣A關(guān)聯(lián)了圖像矢量u和投影矢量p，A=[aij]∈RI×J，其中I為射線總數(shù)，J為圖像的像素總數(shù)，aij為j個像素對第i個投影數(shù)據(jù)的貢獻。CT重建過程是從投影數(shù)據(jù)p中重建圖像u的過程，其實質(zhì)就是線性方程組的求解過程。

由于投影數(shù)據(jù)的不完備性，稀疏角度CT重建是一個不適定問題。對于此類問題，通常引入正則化項，增加先驗條件約束，優(yōu)化模型為：

(2)

1.2 基于TV和AEDS的CT重建算法

基于全變分的CT圖像重建算法。其重建算法的目標(biāo)函數(shù)為：

(3)

(4)

式中，x和y為像素位置；ux,y為圖像在(x,y)處的像素值。該優(yōu)化問題使用梯度下降法進行求解。

AEDS在有限角度的重建中取得了較好的重建效果，該算法是在x和y方向分別引入一個正則項，使圖像趨近于分片常數(shù)的作用，其優(yōu)化模型為：

(5)

(6)

式中，sgn為符號函數(shù)，其屬性為sgn(0)=0。

AEDS經(jīng)驗證，可以有效地保護重建圖像的邊緣，充分利用圖像的稀疏性，彌補以TV作為正則項在重建圖像過程中造成邊緣結(jié)構(gòu)過度平滑，細節(jié)丟失的問題。因此，本文將TV與ADES相結(jié)合，提出一種新的CT重建算法，構(gòu)建目標(biāo)函數(shù)為：

(7)

式中，μ為懲罰參數(shù)；λ1，λ2和λ3分別用來調(diào)節(jié)正則化項在優(yōu)化目標(biāo)函數(shù)中的作用。

1.3 優(yōu)化算法

最優(yōu)化問題(5)可能是非凸的，我們采用一種交替極小化方法進行求解(7)。主要內(nèi)容如下：λ為正常數(shù)，令u(t)表示第t次迭代的結(jié)果，第t+1次的迭代結(jié)果可依次求解以下四個子問題得到：

子問題1：

(8)

子問題2：

(9)

子問題3：

(10)

子問題4：

(11)

子問題1的求解，使用ART進行一次迭代，有：

(12)

子問題2使用梯度下降法進行求解，首先確定梯度下降方向

(13)

(14)

子問題3可以分解為一系列一維問題，首先引入一個輔助變量

g={gij|i=1,2,…,K1,j=1,2,…,K2}

(15)

令gn=gij，其中n=(i-1)K2+j。則子問題3可以轉(zhuǎn)化為

(16)

1)更新重建圖像u，u的求解相當(dāng)于最小化問題

(17)

由于是二次函數(shù)，故有一個全局最小值，得

(18)

2)更新g，其目標(biāo)函數(shù)為

(19)

得g滿足以下條件：

(20)

子問題4的求解同理子問題3。

2 實驗仿真

2.1 實驗設(shè)置

本文以Shepp-Logan Phantom和胸部圖像作為仿真模型，以驗證所提出TV-AEDS算法在稀疏角度CT重建中的性能(圖1)。為了更好地模擬真實的CT場景，在投影正弦圖中加入泊松噪聲，泊松噪聲的量子數(shù)為106，以驗證算法對噪聲的魯棒性。該實驗中X射線源到物體的旋轉(zhuǎn)中心的距離為1000，射線源到探測器的距離為1300，探測器的個數(shù)為392，重建圖像大小為256×256。在稀疏投影角度的重建中，影響重建質(zhì)量的主要因素是投影角度的個數(shù)。因此，在本實驗中，在360°范圍內(nèi)等間隔采樣，采樣間隔分別為12°、6°、4°，得到稀疏投影角度數(shù)目為30、60、90。

圖1 兩種模型，分別為Shepp-Logan Phantom和胸部圖像

在本文的實驗中，迭代算法的參數(shù)以及迭代次數(shù)均通過經(jīng)驗選取，為公平起見，所有參數(shù)都一致。迭代次數(shù)均設(shè)置為1000次，迭代算法中的λ設(shè)置為1。為了進一步定量客觀評價算法的性能，本文使用峰值信噪比(PSNR)[23]、均方根誤差(RMSE)[24]和通用質(zhì)量指標(biāo)(SSIM)[25]來定量地評價重建質(zhì)量。PSNR用來衡量重建圖像質(zhì)量的高低，值越高，重建圖像質(zhì)量越高；反之，重建圖像質(zhì)量越低。RMSE用來衡量重建圖像與參考圖像的像素相似性，值越小，相似性越高；反之，相似性越低。SSIM值越接近1，則重建圖像越接近真實圖像。本文所有算法在Visual Studio 2013和Matlab 2019a上運行[配置Intel(R) Core(TM) i5-10210U CPU @ 1.60GHz 2.11GHz]。

2.2 Shepp-Logan Phantom重建結(jié)果分析

ART-TV算法、ART-AEDS算法、ART-BDTV算法[20]和TV-AEDS算法重建的Shepp-Logan Phantom重建結(jié)果如圖2～圖4所示。由圖可知，稀疏投影角度數(shù)目為30個時，重建圖像最差，圖中存在大量噪聲，且圖像模糊。隨著稀疏投影角度數(shù)目的增加，重建效果越來越好。以稀疏投影角度為60個的重建結(jié)果為例，ART-TV算法和ART-AEDS算法的重建結(jié)果均存在條紋偽影，ART-BDTV算法中有所緩解，效果明顯，而本文提出的算法幾乎完全消除了噪聲，且圖像結(jié)構(gòu)清晰。

圖2 稀疏投影角度數(shù)目為30，分別為ART-TV、ART-AEDS、ART-BDTV和TV-AEDS算法的重建結(jié)果

圖3 稀疏投影角度數(shù)目為60，分別為ART-TV、ART-AEDS、ART-BDTV和TV-AEDS算法的重建結(jié)果

圖4 稀疏投影角度數(shù)目為90，分別為ART-TV、ART-AEDS、ART-BDTV和TV-AEDS算法的重建結(jié)果

四種算法下PSNR、RMSE以及SSIM的定量評價結(jié)果如表1所示。由表可知，本文所提算法的PSNR值最高，RMSE值最低，SSIM值最接近1。因此，本文算法的重建效果相比于其他算法更好。

表1 不同算法在不同稀疏角度下的量化結(jié)果

2.3 胸部圖重建結(jié)果分析

ART算法、ART-TV算法、ART-AEDS算法和TV-AEDS算法重建的胸部圖重建結(jié)果如圖5～圖7所示。幾種算法相比較而言，以60個投影角度為例，ART算法的重建結(jié)果存在大量的噪聲，且圖像較為模糊，結(jié)構(gòu)幾乎分辨不出來。ART-TV算法和ART-AEDS算法的重建結(jié)果噪聲雖然有所減少，但圖像結(jié)構(gòu)仍然很模糊。TV-AEDS算法可以抑制更多的噪聲且重建結(jié)果很清晰，邊緣細節(jié)保護良好，結(jié)構(gòu)完整，可以很好地辨別血管、陰影等情況。

圖5 稀疏投影角度數(shù)目為30，第一行分別為ART、ART-TV、ART-AEDS和TV-AEDS算法的重建結(jié)果；第二行分別為相對應(yīng)的重建圖像方框區(qū)域的放大圖

圖7 稀疏投影角度數(shù)目為90，第一行分別為ART、ART-TV、ART-AEDS和TV-AEDS算法的重建結(jié)果；第二行分別為相對應(yīng)的重建圖像方框區(qū)域的放大圖

四種算法下PSNR、RMSE以及SSIM的定量評價結(jié)果如表2所示。由表可知，本文所提算法的PSNR值最高，RMSE值最低，SSIM值最接近1。因此，本文算法的重建效果相比于其他算法更好。

圖6 稀疏投影角度數(shù)目為60，第一行分別為ART、ART-TV、ART-AEDS和TV-AEDS算法的重建結(jié)果；第二行分別為相對應(yīng)的重建圖像方框區(qū)域的放大圖

表2 不同算法在不同稀疏角度下的量化結(jié)果(胸部圖像)

3 結(jié)論

為了提高稀疏角度CT圖像重建質(zhì)量，本文提出了一種基于TV和AEDS算法的稀疏角度迭代重建算法。通過引入AEDS算法，在水平和垂直兩個方向利用l0范數(shù)的高稀疏性，保留了圖像的更多邊緣結(jié)構(gòu)和細節(jié)，使重建圖像的邊緣以及微小細節(jié)得到了更好的重建效果。實驗結(jié)果表明，與其他算法相比，本文提出的稀疏角度重建算法有較好的性能，不僅可以有效地抑制噪聲，而且能夠很好地保護圖像邊緣信息和細節(jié)特征。