基于全變分與交替保邊擴(kuò)散平滑的稀疏角度CT重建算法

張瑩芳，陳 平，潘晉孝

(1.中北大學(xué) 信息與通信工程學(xué)院，太原 030051；2.中北大學(xué) 信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室，太原 030051)

0 引言

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

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

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

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

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

1 CT重建算法

1.1 CT重建模型

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

Au=p

(1)

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

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

(2)

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

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

(3)

(4)

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

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

(5)

(6)

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

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

(7)

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

1.3 優(yōu)化算法

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

子問(wèn)題1：

(8)

子問(wèn)題2：

(9)

子問(wèn)題3：

(10)

子問(wèn)題4：

(11)

子問(wèn)題1的求解，使用ART進(jìn)行一次迭代，有：

(12)

子問(wèn)題2使用梯度下降法進(jìn)行求解，首先確定梯度下降方向

(13)

(14)

子問(wèn)題3可以分解為一系列一維問(wèn)題，首先引入一個(gè)輔助變量

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

(15)

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

(16)

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

(17)

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

(18)

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

(19)

得g滿足以下條件：

(20)

子問(wèn)題4的求解同理子問(wèn)題3。

2 實(shí)驗(yàn)仿真

2.1 實(shí)驗(yàn)設(shè)置

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

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

在本文的實(shí)驗(yàn)中，迭代算法的參數(shù)以及迭代次數(shù)均通過(guò)經(jīng)驗(yàn)選取，為公平起見(jiàn)，所有參數(shù)都一致。迭代次數(shù)均設(shè)置為1000次，迭代算法中的λ設(shè)置為1。為了進(jìn)一步定量客觀評(píng)價(jià)算法的性能，本文使用峰值信噪比(PSNR)[23]、均方根誤差(RMSE)[24]和通用質(zhì)量指標(biāo)(SSIM)[25]來(lái)定量地評(píng)價(jià)重建質(zhì)量。PSNR用來(lái)衡量重建圖像質(zhì)量的高低，值越高，重建圖像質(zhì)量越高；反之，重建圖像質(zhì)量越低。RMSE用來(lái)衡量重建圖像與參考圖像的像素相似性，值越小，相似性越高；反之，相似性越低。SSIM值越接近1，則重建圖像越接近真實(shí)圖像。本文所有算法在Visual Studio 2013和Matlab 2019a上運(yùn)行[配置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個(gè)時(shí)，重建圖像最差，圖中存在大量噪聲，且圖像模糊。隨著稀疏投影角度數(shù)目的增加，重建效果越來(lái)越好。以稀疏投影角度為60個(gè)的重建結(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的定量評(píng)價(jià)結(jié)果如表1所示。由表可知，本文所提算法的PSNR值最高，RMSE值最低，SSIM值最接近1。因此，本文算法的重建效果相比于其他算法更好。

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

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

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

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

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

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

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

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

3 結(jié)論

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