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

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

0 引 言

計(jì)算機(jī)斷層掃描技術(shù)(Computed Tomography，CT)作為20世紀(jì)后期影像學(xué)的杰出代表，早已廣泛應(yīng)用于健康普查中[1]. 由于臨床中常規(guī)劑量的CT掃描具有潛在的致癌風(fēng)險(xiǎn)，因此，“低劑量”已經(jīng)成為CT研究的主流方向. 近年來，一些顯著降低輻射劑量的研究主要聚焦于稀疏投影數(shù)據(jù)重建方面. 但稀疏重建由于投影角度的減少，導(dǎo)致了重建圖像往往存在明顯的條形偽影.

目前，解決該問題的一類主要方法是引入全變差(Total Variation，TV)稀疏性正則項(xiàng)作為約束條件[2-3]. 在此基礎(chǔ)上，相關(guān)改進(jìn)算法也得到了深入研究. 2009年Sidky[4]提出利用稀疏性更強(qiáng)的Lp范數(shù)來代替TV中的L1范數(shù)，獲得了更加精確的重建圖像； 針對(duì)傳統(tǒng)TV正則項(xiàng)缺乏方向信息，易模糊低對(duì)比度邊緣，2010年Jin等人[5]提出各向異性全變差(Anisotropic TV，ATV)，同時(shí)利用圖像稀疏性和邊緣方向信息，在一定程度上來保護(hù)邊緣； Bayram[6]等人于2012年提出一種方向全變差(Directional TV，DTV)圖像去噪模型，通過增加某一主導(dǎo)方向的權(quán)重，顯著增強(qiáng)該方向的邊緣. 由于權(quán)重及方向參數(shù)是固定的，與ATV一樣，該方法丟失了大量其他方向上的邊緣信息； 為了克服該缺陷，Tao等人[7]利用估計(jì)的結(jié)構(gòu)張量場計(jì)算圖像方向結(jié)構(gòu)，引入了多方向的邊緣信息，進(jìn)一步提高了重建質(zhì)量.

鑒于此，本文提出了一種基于結(jié)構(gòu)張量的自適應(yīng)方向全變差稀疏角度CT重建算法，以下簡稱SADTVp.

1 背景知識(shí)

1.1 基于TV的CT圖像重建模型

稀疏角度CT重建問題是一個(gè)典型的病態(tài)逆問題，常用圖像f的全變差作為其稀疏正則項(xiàng)，其優(yōu)化模型為

(1)

(2)

式中： Δ1和Δ2分別表示水平和垂直離散算子.

1.2 DTV正則項(xiàng)

將式(2)改寫為

(3)

式中：B2是L2范數(shù)的單位球，Δf=(Δ1f,Δ2f)T.

傳統(tǒng)TV模型中Δf的各分量表現(xiàn)為各向同性，缺乏方向信息，這將降低對(duì)低對(duì)比度邊緣的檢測能力. 為了克服這一不足，Bayram和Kamasak利用與圖像梯度分布相對(duì)應(yīng)的橢圓Eα,θ替換B2，提出了如下DTV正則項(xiàng)模型[6]

(4)

式中：Eα,θ=RθΛαB2，Rθ、Λα分別表示旋轉(zhuǎn)矩陣與權(quán)重矩陣，用來刻畫圖像中特征結(jié)構(gòu)的方向及該方向的權(quán)重[8]，表達(dá)式為

(5)

(6)

式中：RTθ=R-θ，θ為邊緣方向與x軸夾角，α>1，表示沿邊緣方向的尺度參數(shù).

DTV模型中θ與α參數(shù)的選取是固定的，只對(duì)沿角度θ方向的邊緣信息敏感，這將丟失大量其他方向的邊緣信息.

2 本文算法

2.1 算法思想

為了能夠充分利用各個(gè)方向的邊緣信息，受文獻(xiàn)[7,9]啟發(fā)，本文引入結(jié)構(gòu)張量，在重建過程中自適應(yīng)地獲取θ與α值.

結(jié)構(gòu)張量常用來刻畫圖像局部區(qū)域的結(jié)構(gòu)方向信息，與圖像梯度相比，提取的邊緣信息更加精確，其定義為

(7)

式中：uδ表示對(duì)圖像進(jìn)行參數(shù)為δ的高斯卷積，本文采用雙邊濾波代替高斯平滑，力求達(dá)到去噪與保留邊緣的平衡；gρ為尺度為ρ的高斯核，引入局部鄰域信息.

由于Jρ(u)是對(duì)稱且半正定的二維矩陣，對(duì)圖像邊緣的方向估計(jì)和結(jié)構(gòu)分析等可以通過對(duì)其進(jìn)行特征分解實(shí)現(xiàn)[9]. 兩個(gè)非負(fù)特征值為

(8)

對(duì)應(yīng)的特征向量為

(9)

分析式(5)可知，當(dāng)a?1時(shí)，等價(jià)于僅沿角度為θ的邊緣方向進(jìn)行正則化，故本文將權(quán)重矩陣改寫為

(10)

其中，0≤α′≤1.

基于結(jié)構(gòu)張量的特點(diǎn)，本文將圖像中每一像素點(diǎn)的θ和α′分別定義為

(11)

當(dāng)處于邊緣區(qū)域時(shí)α′≈0，圖像僅沿邊緣方向進(jìn)行各向異性擴(kuò)散； 處于平坦區(qū)域時(shí)α′≈1，圖像進(jìn)行各向同性擴(kuò)散. 然而由于重建過程中的中間圖像存在偽影及噪聲，直接利用結(jié)構(gòu)張量得到的初始邊緣圖通常較粗糙，不適合用來約束DTV正則化. 本文對(duì)初始邊緣圖進(jìn)行邊緣細(xì)化處理，剔除偽邊緣點(diǎn)后，再更新式(11)，即得到較準(zhǔn)確的自適應(yīng)參數(shù).

2.2 算法模型

綜合計(jì)算式(4)與式(6)～(11)，同時(shí)考慮到，Lp范數(shù)(0

(12)

式中：f=[f1,f2,…,fN]T為圖像矢量.

因此，重建模型為

(13)

對(duì)應(yīng)的離散形式為

(14)

3 模型求解

鑒于式(14)中的SADTVp范數(shù)正則項(xiàng)是非凸且非光滑的，本文基于宋潔等人提出的非凸模型求解算法[11]，采用交替方向乘子法ADMM結(jié)合廣義軟閾值算法進(jìn)行求解. 引入輔助變量wi=Λα′aiR-θaiΔif, 可以將式(14)轉(zhuǎn)化為約束優(yōu)化問題，即

s.t.wi=Λα′aiR-θaiΔif,

(15)

其中，wi∈R2，i=1,2,…,N. 為了便于計(jì)算，令Λα′aiR-θaiΔi=Ti，則式(15)的增廣拉格朗日函數(shù)可表示為

(16)

式中：λ為拉格朗日乘子；μ為懲罰參數(shù).

利用ADMM最小化式(16)相當(dāng)于在固定一個(gè)參數(shù)不變的情況下交替求解w與f兩個(gè)子問題并更新λ的值[12-14].

3.1 w子問題求解

已知λ，保持f不變，最小化w可以表示為

(17)

式中：k表示迭代次數(shù). 為解決上述Lp范數(shù)最小化問題，文獻(xiàn)[15]提出一種廣義軟閾值函數(shù)，

TGSTP(w,λ)=

(18)

(19)

3.2 f子問題求解

已知λ值，保持w不變，f子問題求解可以表示為

(20)

顯然，最小化f是一個(gè)最小二乘問題，此處利用梯度下降法求解. 首先計(jì)算目標(biāo)函數(shù)的偏導(dǎo)數(shù)為

dk=β[AT(Afk-p)]+

(21)

從而得到

fk+1=fk-αdk,

(22)

式中：α為下降步長.

3.3 算法描述

綜上所述，SADTVp算法步驟概括如下：

1) 輸入初始圖像f0，初始化參數(shù)λ0,w0,Rθ0,Λa′0,00,μ>0,n,ε.

2) 利用式(19)更新w.

3) 利用式(22)更新f.

4) 更新λ=λ-μ(w-Tf).

5) 利用式(6)～(11)計(jì)算Rθ和Λα′，得到初始邊緣圖.

6) 對(duì)初始邊緣圖進(jìn)行邊緣細(xì)化.

7) 利用細(xì)化后的邊緣圖更新Rθ和Λα′.

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

本文采用圖1(a) 所示的大小為256×256的Shepp-Logan頭模作為實(shí)驗(yàn)對(duì)象. 在0°～180°范圍均勻采集60個(gè)角度，探測器個(gè)數(shù)為260個(gè). 實(shí)驗(yàn)中，在濾波反投影(Filtered Back Projection，F(xiàn)BP)重建過程中，利用零投影數(shù)據(jù)將空氣區(qū)域修正為零，將非局部均值濾波后的重建圖像作為迭代的初始圖像f0，如圖1(b) 所示.

圖1 Shepp-Logan模型及初始圖Fig.1 Shepp-Logan model and the initial image

為驗(yàn)證算法的有效性，選用FBP、 TV、 ATV[5]、 各向異性邊緣檢測引導(dǎo)的全變差CT重建算法EGTV[16]以及文獻(xiàn)[11]中基于Lp范數(shù)的TV算法作為對(duì)比算法，在無噪聲的情況下進(jìn)行算法驗(yàn)證. 同時(shí)，將所提算法的Lp范數(shù)在p=1時(shí)定義為SADTV，參與對(duì)比分析.

EGTV算法依據(jù)所在文獻(xiàn)，采用凸集投影法求解，其余算法均采用ADMM方法求解. 經(jīng)試驗(yàn)，相關(guān)參數(shù)設(shè)置如下：

μ=26，β=0.6，ε=1×10-3，λ0=0，w0=0, 對(duì)比算法中的其他參數(shù)均參考其所在文獻(xiàn).Lp范數(shù)中p=0.7.

此外，為進(jìn)一步評(píng)價(jià)所提算法的重建性能，采用歸一化平均絕對(duì)距離判據(jù)NAAD、 結(jié)構(gòu)相似度SSIM來定量評(píng)價(jià)各算法重建效果.

(23)

(24)

式中：ti,j和ri,j分別表示原始圖像和重建后圖像中第i行、j列的像素密度，μr、μt、δr、δt、δr,t分別表示重建后圖像與原始圖像的密度平均值、 標(biāo)準(zhǔn)差及協(xié)方差.

圖2 為TV算法在60個(gè)采樣角度下重建圖像的NAAD、SSIM與迭代次數(shù)關(guān)系曲線圖. 可以看出，迭代次數(shù)為70次以后曲線趨于平穩(wěn)，因此迭代次數(shù)設(shè)定為70次. 圖3 為不同算法在60個(gè)角度下迭代70次的重建結(jié)果. 同時(shí)，選取如圖3(a)中矩形框標(biāo)注的局部區(qū)域放大顯示，如圖4 所示.

圖2 TV算法在60個(gè)角度下重建圖像的NAAD、 SSIM與迭代次數(shù)關(guān)系曲線Fig.2 Iterative curves of NAAD and SSIM of images reconstructed by TV from 60 sparse-view

圖3 不同算法在60個(gè)角度下迭代70次的重建圖像Fig.3 Reconstructed images based on 60 angles after 70 times of iteration using different algorithms

圖4 矩形標(biāo)注區(qū)域的放大圖Fig.4 Enlarged drawing of the rectangular label area

觀察分析圖3 及圖4 可以發(fā)現(xiàn)，F(xiàn)BP算法重建圖像含有大量條形偽影，邊緣細(xì)節(jié)模糊； TV算法重建效果明顯優(yōu)于FBP算法，但在垂直方向含有少量偽影，同樣存在邊緣細(xì)節(jié)模糊的現(xiàn)象； ATV算法重建圖像則呈現(xiàn)出垂直、 水平方向偽影； 與TV算法相比，文獻(xiàn)[11]中算法、 SADTV算法和EGTV算法的重建質(zhì)量均有所提升. 文獻(xiàn)[11]中算法雖較好抑制了偽影，但相鄰邊緣仍有連接； SADTV算法和EGTV算法的重建圖像的邊緣清晰度有所提升，但垂直方向仍存在偽影. 相比而言，本文算法重建圖像具有較好的一致性，邊緣細(xì)節(jié)清晰，與原始圖像的接近程度較高.

圖5 為原圖及六種迭代重建算法圖像在水平、 垂直方向的剖面輪廓圖.

圖5 重建圖像剖面輪廓對(duì)比Fig.5 Profile contour contrast of reconstructed images

由圖5 可以看出，本文算法與原始圖像輪廓最為吻合，尤其在邊緣跳躍部分，改善了其他算法邊緣鈍化現(xiàn)象. 各算法迭代70次的質(zhì)量參數(shù)見表1. 綜合來看，七種算法中本文算法的SSIM最高，NAAD最小，算法的有效性得到驗(yàn)證.

表1 七種重建算法的質(zhì)量參數(shù)對(duì)比

5 結(jié) 論

針對(duì)稀疏角度CT重建中全變差模型的局限性，提出將具有方向先驗(yàn)的方向全變差與更具稀疏性的Lp范數(shù)結(jié)合. 基于結(jié)構(gòu)張量提取局部結(jié)構(gòu)信息的優(yōu)勢，自適應(yīng)地提取多個(gè)方向的邊緣信息. 實(shí)驗(yàn)表明，與FBP算法、 TV算法、 ATV算法、 文獻(xiàn)[11]中算法、 SADTV算法、 EDTV算法相比，所提SADTVp算法能有效抑制條形偽影且圖像邊緣輪廓完整清晰，可以提高稀疏重建性能. SADTVp算法重建效果依賴于結(jié)構(gòu)張量邊緣檢測的精度，本文下一步將研究更優(yōu)的邊緣檢測方法，提取更加精細(xì)的邊緣，以進(jìn)一步提升算法性能.