超拉普拉斯重疊組稀疏先驗(yàn)的稀疏角度CT 重建

齊子文，孔慧華*，李佳欣，潘晉孝

1 中北大學(xué)數(shù)學(xué)學(xué)院，山西 太原 030051；

2 信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051

1 引言

X 射線計(jì)算機(jī)斷層掃描(computed tomography,CT)成像技術(shù)可以獲取人體和物體的內(nèi)部結(jié)構(gòu)信息，已廣泛應(yīng)用于工業(yè)無損檢測(cè)、醫(yī)學(xué)臨床診斷等領(lǐng)域[1]。然而，在臨床檢查時(shí)額外的X 射線輻射可能導(dǎo)致癌癥和其他遺傳病變[2]。因此，降低輻射劑量和重建高質(zhì)量圖像是目前的研究重點(diǎn)。降低X 射線輻射量的方法主要分兩方面：第一種方法是減少每個(gè)投影中的X 射線曝光，第二種方法是減少投影采樣數(shù)，即稀疏投影[3]。第一種方法會(huì)產(chǎn)生含有噪聲的投影，重建圖像含有大量的噪聲和偽影[4]。第二種方法是應(yīng)用稀疏角度CT，但由于投影數(shù)據(jù)的采樣不足，重建圖像會(huì)產(chǎn)生嚴(yán)重的偽影[5]。

一些傳統(tǒng)的CT 圖像重建算法比如濾波反投影(filtered back projection,FBP)算法和聯(lián)合代數(shù)重建(simultaneous algebraic reconstruction technique,SART)算法，這兩個(gè)基礎(chǔ)算法無法對(duì)不完整的投影數(shù)據(jù)重建出高質(zhì)量的圖像。壓縮感知(Compressed sensing,CS)理論的出現(xiàn)使得圖像處理中的不適定問題可以有效地解決[6]?；贑S 的圖像重建算法主要利用圖像的先驗(yàn)知識(shí)，特別是圖像的稀疏性，如全變分[7]、字典學(xué)習(xí)[8]、小波變換等[9]。近年來，全變分(total variational,TV)算法被廣泛應(yīng)用于不同場(chǎng)景下的CT 重建，該算法可以很好地平滑圖像以達(dá)到去噪的目的[10]。Sidky 等將TV 最小化引入到不完全投影數(shù)據(jù)重建中，提出了一種基于投影到凸集(projections onto convex sets,POCS)的高效迭代重建算法TV-POCS[11]。雖然TV 正則項(xiàng)能夠很好地恢復(fù)圖像的邊緣信息，但是過渡的平滑圖像會(huì)導(dǎo)致圖像產(chǎn)生階梯偽影。為了克服這一現(xiàn)象，Zhang 等將非零像素稀疏性和圖像連續(xù)性假設(shè)兩個(gè)正則項(xiàng)來代替TV 正則項(xiàng)，提出了一種基于大規(guī)模線性規(guī)劃的稀疏角度層析重建算法[12]。連祥媛等提出一種多通道聯(lián)合的廣義全變分 (Total Generalized Variational,TGV)的能譜CT 迭代重建算法，利用通道之間的相關(guān)性來消除重建過程中產(chǎn)生的噪聲和階梯偽影[13]。Peng 等將組稀疏(group sparse representation,GSR)正則化引入稀疏角度CT 重建，利用圖像中相似補(bǔ)丁構(gòu)成的圖像組來表示GSR 的基本單元，并且成功將其應(yīng)用于稀疏角度CT 重建中[14]。Jon 等將重疊組稀疏(overlapping group sparsity,OGS)和超拉普拉斯先驗(yàn)(hyper-Laplacian,HL)相結(jié)合，提出了一種基于重疊組稀疏超拉普拉斯先驗(yàn)(overlapping group sparsity on hyper-Laplacian,OGSHL)的圖像去噪算法[15]。

本文將OGS 正則項(xiàng)和HL 正則項(xiàng)應(yīng)用于稀疏角度CT 圖像重建。相比于傳統(tǒng)的TV 正則項(xiàng)，OGS 正則項(xiàng)改進(jìn)圖像梯度的稀疏性，由于重疊組采用圖像梯度的結(jié)構(gòu)信息作為圖像梯度稀疏性的度量標(biāo)準(zhǔn)，因此該算法能夠在恢復(fù)圖像邊緣的同時(shí)消除階梯偽影，基于HL 先驗(yàn)可以很好地近似圖像梯度的重尾分布，對(duì)于重建圖像細(xì)節(jié)的保留有著更好的效果。

2 理論推導(dǎo)

2.1 基于TV 的稀疏角度CT 重建算法

稀疏角度CT 圖像重建與傳統(tǒng)CT 圖像重建的過程基本相同，設(shè)稀疏角度CT 重建中需要重建的圖像為f，從探測(cè)器上獲得的投影數(shù)據(jù)為p，數(shù)據(jù)采集的投影矩陣用A表示，稀疏角度CT 重建問題可表示為

在稀疏角度CT 重建的情況下，由于投影采樣不足會(huì)使重建圖像產(chǎn)生噪聲和偽影。為了重建高質(zhì)量圖像，使用CS 理論將圖像進(jìn)行稀疏變換，例如對(duì)重建圖像進(jìn)行離散梯度變換(discrete gradient transform,DGT)，則任意圖像的DGT 定義為

則對(duì)于任意圖像f，其TV 表達(dá)式為

其中?x，?y分別表示水平方向和豎直方向上的梯度算子，其表達(dá)式分別為

故基于TV 正則項(xiàng)的稀疏角度CT 重建模型為

其中：第一項(xiàng)為數(shù)據(jù)保真項(xiàng)，A為投影矩陣，p為投影數(shù)據(jù)，μ為保真項(xiàng)系數(shù)。第二項(xiàng)為TV 正則項(xiàng)，λ為正則項(xiàng)系數(shù)。

2.2 基于OGS-TV 的稀疏角度CT 重建算法

TV 正則項(xiàng)已被廣泛應(yīng)用，其雖然能夠保留圖像的銳利邊緣，但是過渡的平滑圖像會(huì)導(dǎo)致圖像產(chǎn)生階梯偽影。為了進(jìn)一步提升圖像的質(zhì)量，消除重建中產(chǎn)生的階梯偽影，Liu 等將一維信號(hào)的OGS 推廣到二維圖像中并成功應(yīng)用于圖像去噪[16]。

為了描述圖像梯度的結(jié)構(gòu)稀疏性，首先定義一個(gè)二維圖像的像素組，形式如

其中：d(i,j),K是矩陣中所有元素組成的列向量，即：d(i,j),K=。根據(jù)式(8)可知當(dāng)K=1 時(shí)，該正則項(xiàng)就轉(zhuǎn)變?yōu)閭鹘y(tǒng)的TV 正則項(xiàng)，則OGS-TV 正則項(xiàng)可定義為

由于OGS 正則項(xiàng)考慮了每個(gè)像素點(diǎn)梯度的K×K方形領(lǐng)域內(nèi)的信息，是一個(gè)非局部的稀疏先驗(yàn)，因此OGS-TV 正則項(xiàng)不僅擁有和TV 正則項(xiàng)一樣的圖像邊緣保留能力，還能克服TV 重建中圖像邊緣產(chǎn)生的階梯偽影。故基于OGS-TV 的稀疏角度CT 重建算法為

其中：第一項(xiàng)為數(shù)據(jù)保真項(xiàng)，第二項(xiàng)為OGS-TV 正則項(xiàng)，λ為正則項(xiàng)系數(shù)。

2.3 基于OGS-HL 的稀疏角度CT 重建算法

OGS-TV 雖然可以克服圖像的階梯偽影，但它難以恢復(fù)圖像中較為復(fù)雜的紋理和細(xì)節(jié)。為進(jìn)一步擴(kuò)展OGS-TV，Kong 等強(qiáng)調(diào)圖像的梯度服從重尾分布[17]，并且超拉普拉斯算子先驗(yàn)將會(huì)使得重尾分布獲得更好的近似，提出的HL 先驗(yàn)的模型為

其中：‖·‖q表示擬范數(shù)lq，0＜q＜1。q為超拉普拉斯參數(shù)，通常參數(shù)q的范圍為 0.5≤q≤1。根據(jù)式(11)，，將HL 先驗(yàn)代入OGS正則項(xiàng)得

根據(jù)式(12)，當(dāng)超拉普拉斯參數(shù)q=1時(shí)，該正則項(xiàng)就轉(zhuǎn)變?yōu)镺GS-TV 正則項(xiàng)。OGS-HL 的正則項(xiàng)可表示為

故基于OGS-HL 的稀疏角度CT 重建算法，算法模型為

其中：第一項(xiàng)為數(shù)據(jù)保真項(xiàng)，第二項(xiàng)為OGS-HL 正則項(xiàng)，λ為正則項(xiàng)系數(shù)。

3 模型求解

OGS-HL 算法的模型是一個(gè)綜合算法模型，TV，OGS-TV，OGS-HL 算法可以通過對(duì)模型設(shè)置不同的參數(shù)值來分別表示，如表1 所示。

表1 不同算法的參數(shù)設(shè)置Table 1 Parameter settings for different algorithms

由于算法模型包含的矩陣數(shù)據(jù)量大并且復(fù)雜，因此本文采用交替方向乘子 (alternating direction method of multipliers,ADMM)算法[18]，將模型轉(zhuǎn)化為多個(gè)子問題來求解，式(14)引入變量u，將模型轉(zhuǎn)化為

上述優(yōu)化函數(shù)可分解為2 個(gè)子問題來分別更新f和u：

其中：式(17)采用梯度下降法來更新f：

式(18)進(jìn)一步采用ADMM 算法來求解，引入輔助變量x1和x2，引入拉格朗日乘子w1，w2和懲罰參數(shù)δ，轉(zhuǎn)化為以下包含多個(gè)變量的增廣拉格朗日函數(shù)來更新：

Step1：更新x1和x2，根據(jù)式(20)，x1和x2的子問題可以由式(21)和式(22)來更新：

該子問題的結(jié)構(gòu)比較復(fù)雜，故采用主分量最小化法(majorization minimization,MM)來求解上述最小化問題[19]。由于2 個(gè)子問題是相同的結(jié)構(gòu)，下面只闡述式(21)的更新步驟。

為了方便后續(xù)計(jì)算，式(21)可簡(jiǎn)化為

首先，需找到φOH(x1)的優(yōu)化目標(biāo)函數(shù)，根據(jù)均值不等式，可以得到：

其中：a(i,j),K是 φOH(x1)中的任意一項(xiàng)，對(duì)于任意的a(i,j),K，x1(i,j),K≠0，并且a(i,j),K=x1(i,j),K時(shí)，對(duì)每一個(gè)稀疏組，代入式(24)可以得到 φOH(x1)的優(yōu)化目標(biāo)函數(shù)：

根據(jù)式 (24)能夠證明 ψ(x1,a)≥φ(x1)，并且ψ(a,a)=φ(a)。通過簡(jiǎn)單的計(jì)算，ψ(x1,a)可以寫為

其中：C 是一個(gè)常數(shù)，Λ(a)是對(duì)角矩陣，其對(duì)角的元素為

其中：r,t=1,2,...,n,l=1,2,...,n2，Λ(a)是通過二維卷積運(yùn)算來求解，則式(23)中R(x1)的優(yōu)化目標(biāo)函數(shù)如：

式(28)可以得出 ξ(x1,a)≥R(x1)，ξ(a,a)=R(a)。故 ξ(x1,a)也滿足MM 算法的前提條件，為了極小化R(x1)，使用MM 算法迭代，初始化，反復(fù)極小化替代函數(shù)，得到

根據(jù)歐拉-拉格朗日方程，式(29)可寫為

其中：E是一個(gè)全1 的單位矩陣，S(x1)=diag(|x1|2q-2)。

Step2：更新w1和w2，拉格朗日乘子w1和w2的更新公式為

Step3：更新u，u的子問題的更新公式為

上述u子問題的優(yōu)化目標(biāo)是二次泛函，故可以求解下述等效正態(tài)方程：

顯然，式(34)含有二維卷積運(yùn)算，因此不能直接求出u子問題的最優(yōu)解，空域中的卷積運(yùn)算在頻域中會(huì)轉(zhuǎn)化為普通的乘法運(yùn)算，故采用二維快速傅里葉變換與逆變換來求解：

其中：F和F-1分別是快速傅里葉變換和快速傅里葉逆變換。

基于OGS-HL 的稀疏角度CT 重建算法步驟：

4 實(shí)驗(yàn)結(jié)果

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

為驗(yàn)證OGS-HL 算法在稀疏角度CT 重建下的效果，本文選取其他經(jīng)典的算法與OGS-HL 算法作對(duì)比。實(shí)驗(yàn)采用仿真模型和真實(shí)模型來驗(yàn)證算法。實(shí)驗(yàn)配置為：11th Gen Intel (R) Core (TM) i5-11260H @ 2.60 GHz 的CPU，Nvidia RTX3050ti (4 GB)的GPU，16 GB 的內(nèi)存，實(shí)驗(yàn)軟件為MATLAB R2019b。本文算法與選取的對(duì)比算法均在MATLAB 和使用MEX 函數(shù)編譯的C++的混合模式中實(shí)現(xiàn)，接口在MATLAB中實(shí)現(xiàn)。

4.2 仿真小鼠實(shí)驗(yàn)

仿真實(shí)驗(yàn)采用MOBY 軟件生成的小鼠胸腔截面作為測(cè)試模型，在血液中注入1.2%的碘造影劑，胸腔模型尺寸為20 mm×20 mm，分辨率為512×512，如圖1 所示。

圖1 仿真小鼠胸腔模型和碘造影劑Fig.1 The simulated mouse thorax phantom and the iodine contrast agent

投影數(shù)據(jù)由等距扇形束幾何采集，X 射線源到旋轉(zhuǎn)中心的距離為100 mm，物體半徑為10 mm，每個(gè)探測(cè)器長度為20 mm，探測(cè)器有320 個(gè)單元，每個(gè)單元的長度為0.0625 mm。每條X 射線發(fā)射的光子數(shù)為5×104，計(jì)算沿每條X 射線路徑的預(yù)期光子數(shù)，為了模擬數(shù)據(jù)噪聲，根據(jù)泊松分布生成隨機(jī)數(shù)，其中方差為上述期望的光子數(shù)，然后進(jìn)行對(duì)數(shù)運(yùn)算，得到無噪聲和有噪聲的投影數(shù)據(jù)。使用分裂布雷格曼(Split-Bregman,SB)算法的重建的無噪聲圖像作為參考圖像[20]。在360°范圍內(nèi)進(jìn)行等間隔采樣，設(shè)置的采樣間隔角度分別為6°、4°、3°和2°，得到的稀疏投影角度個(gè)數(shù)為60、90、120 和180。本節(jié)實(shí)驗(yàn)采用歸一化均方根誤差(NRMSE)、峰值信噪比(PSNR)和結(jié)構(gòu)相似度指數(shù)(SSIM)來評(píng)價(jià)噪聲投影數(shù)據(jù)中重建圖像的質(zhì)量。

在OGS-TV 算法中，經(jīng)過多次實(shí)驗(yàn)，參數(shù)為μ=50，λ=2×10-4，δ=1，q=1，K=3。在OGS-HL算法中，參數(shù)為 μ=50，λ=2×10-4，δ=1，q=0.8，K=3。

圖2 展示了6 種不同算法下不同角度的仿真小鼠重建圖像，圖3 展現(xiàn)了仿真小鼠骨的部分ROI 區(qū)域。從圖2 可以看出，F(xiàn)BP 和SART 兩種算法在四個(gè)不同的角度下重建的圖像中含有大量的偽影和噪聲。從圖3 的ROI 區(qū)域中仿真小鼠的骨切面可以看出TV 算法重建的圖像出現(xiàn)階梯偽影，使得骨切面的邊緣變得不清晰。TGV 算法雖然能夠克服TV 算法出現(xiàn)的問題，但對(duì)60、90 和120 角度下骨切面周圍產(chǎn)生的條紋偽影去除效果較差。OGS-TV 算法在稀疏角度下重建的圖像都有較大的改進(jìn)，但對(duì)噪聲的抑制效果一般。OGS-HL 算法不僅能夠消除稀疏角度下產(chǎn)生的偽影，并且能夠很好地抑制噪聲，對(duì)小鼠骨切面周圍的噪聲消除以及骨切面的邊緣保護(hù)有著良好的效果。

圖2 仿真小鼠重建結(jié)果: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。從上到下依次為投影幅數(shù)為60、90、120 和180 的重建圖像Fig.2 The reconstruction results of the simulated mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the reconstructed images are from 60,90,120 and 180

圖3 仿真小鼠重建圖像ROI 區(qū)域: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。從上到下依次為投影幅數(shù)為60、90、120 和180 的ROI 區(qū)域Fig.3 The reconstruction results of the simulated mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the ROI regions are 60,90,120 and 180

圖4 繪制了五種算法在重建圖像過程中在四個(gè)不同角度下NRMSE、PSNR 和SSIM 評(píng)價(jià)指標(biāo)的收斂曲線圖。從圖4 可以看出，OGS-HL 算法對(duì)于圖像重建的效果對(duì)比其他算法有所提升，其在180 角度下的提升較小，在120 角度下重建的效果提升，在60 和90 角度下的重建效果提升較為明顯。仿真試驗(yàn)結(jié)果表明，OGS-HL 算法在稀疏角度重建的條件下會(huì)保持良好的優(yōu)勢(shì)，能夠提升重建圖像的質(zhì)量。

4.3 臨床小鼠實(shí)驗(yàn)

為了進(jìn)一步說明OGS-HL 算法的有效性，本節(jié)實(shí)驗(yàn)使用MARS (Medipix All Resolution System)微型CT 上采集的來自真實(shí)臨床前小鼠的投影數(shù)據(jù)。重建的CT 圖像覆蓋面積為18.41 mm×18.41 mm，分辨率為512×512。投影數(shù)據(jù)采用等距扇形束幾何采集，X射線源到旋轉(zhuǎn)中心的距離為158 mm，到探測(cè)器的距離為255 mm。在360°范圍內(nèi)進(jìn)行等間隔采樣，設(shè)置的采樣間隔角度分別為6°、4°、3°和2°，得到的稀疏投影角度個(gè)數(shù)為60、90、120 和180。

在OGS-TV 算法中，經(jīng)過多次實(shí)驗(yàn)，參數(shù)為μ=50，λ=2×10-4δ=0.1，q=1，K=3。在OGS-HL算法中，參數(shù)為 μ=50，λ=2×10-4，δ=0.1，q=0.8，K=3。

圖5 展示了6 種不同算法下不同角度的臨床真實(shí)小鼠重建圖像，圖6 展現(xiàn)了臨床真實(shí)小鼠骨的部分ROI 區(qū)域。從圖5 可以看出，每個(gè)算法重建的效果與以上仿真小鼠重建實(shí)驗(yàn)的效果基本相同，在四個(gè)不同的角度下FBP 算法和SART 算法重建過程產(chǎn)生了大量的偽影和噪聲。從圖6 的ROI 區(qū)域看出，TV 算法下小鼠骨切面的邊緣出現(xiàn)偽影，且含有少量的噪聲，TGV 算法在60、90 角度和120 角度下的噪聲去除并不明顯，OGS-TV 算法也含有少量的噪聲，OGS-HL對(duì)于骨切面周圍噪聲的消除和邊緣的保護(hù)有良好的效果。從實(shí)驗(yàn)結(jié)果和ROI 區(qū)域可以看出，OGS-HL 算法在稀疏角度條件下的重建效果優(yōu)于其他算法。

圖5 臨床小鼠模型重建結(jié)果: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。從上到下依次為投影幅數(shù)為60、90、120 和180 的重建圖像Fig.5 The reconstruction results of the clinical mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the reconstructed images are from 60,90,120 and 180

圖6 臨床小鼠模型重建圖像ROI 區(qū)域: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。從上到下依次為投影幅數(shù)為60、90、120 和180 的ROI 區(qū)域Fig.6 The reconstruction results of the clinical mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the ROI regions are 60,90,120,and 180

5 結(jié)論

為了進(jìn)一步提升稀疏角度CT 圖像的重建質(zhì)量，本文將OGS 和HL 應(yīng)用于稀疏角度CT 圖重建，得到OGS-HL 稀疏角度CT 迭代重建算法。相比于傳統(tǒng)的TV 正則項(xiàng)，OGS 正則項(xiàng)考慮每個(gè)像素點(diǎn)梯度重疊組的信息，度量圖像梯度稀疏性，克服了圖像重建過程中產(chǎn)生的階梯偽影，取得了不錯(cuò)的重建效果。圖像的梯度服從重尾分布，HL 先驗(yàn)可以很好地近似圖像梯度的重尾分布，達(dá)到抑制噪聲的效果。通過以上仿真小鼠和臨床小鼠兩個(gè)實(shí)驗(yàn)結(jié)果表明，OGS-HL 算法對(duì)于稀疏角度CT 圖像重建的圖像質(zhì)量提升明顯，具有較好的魯棒性。本文是基于模型的稀疏CT 重建方法，側(cè)重于物理模型和數(shù)學(xué)算法，而基于深度學(xué)習(xí)的稀疏CT 重建方法則注重于數(shù)據(jù)驅(qū)動(dòng)和自主學(xué)習(xí)能力，更側(cè)重于大量的CT 數(shù)據(jù)作為訓(xùn)練樣本，通過學(xué)習(xí)數(shù)據(jù)中的信息來實(shí)現(xiàn)重建。由于本文進(jìn)行的實(shí)驗(yàn)沒有大量的訓(xùn)練數(shù)據(jù)集，故采用基于模型的算法來實(shí)現(xiàn)重建。在參數(shù)選擇上，正則項(xiàng)參數(shù) λ和懲罰參數(shù) δ的選擇對(duì)重建圖像的效果極為重要，但是，本文這兩個(gè)參數(shù)的選取為經(jīng)驗(yàn)選取，導(dǎo)致進(jìn)行實(shí)驗(yàn)時(shí)耗費(fèi)的時(shí)間較長，之后的研究會(huì)進(jìn)一步改進(jìn)模型，使算法模型更加靈活簡(jiǎn)便。