張 萌，梁亞星，陳 燕，劉佳鑫，桂志國,，張 權,

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

0 引 言

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

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

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

1 背景知識

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

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

(1)

(2)

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

1.2 DTV正則項

將式(2)改寫為

(3)

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

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

(4)

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

(5)

(6)

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

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

2 本文算法

2.1 算法思想

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

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

(7)

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

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

(8)

對應的特征向量為

(9)

分析式(5)可知，當a?1時，等價于僅沿角度為θ的邊緣方向進行正則化，故本文將權重矩陣改寫為

(10)

其中，0≤α′≤1.

基于結構張量的特點，本文將圖像中每一像素點的θ和α′分別定義為

(11)

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

2.2 算法模型

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

(12)

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

因此，重建模型為

(13)

對應的離散形式為

(14)

3 模型求解

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

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

(15)

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

(16)

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

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

3.1 w子問題求解

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

(17)

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

TGSTP(w,λ)=

(18)

(19)

3.2 f子問題求解

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

(20)

顯然，最小化f是一個最小二乘問題，此處利用梯度下降法求解. 首先計算目標函數(shù)的偏導數(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)計算Rθ和Λα′，得到初始邊緣圖.

6) 對初始邊緣圖進行邊緣細化.

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

4 實驗結果與分析

本文采用圖1(a) 所示的大小為256×256的Shepp-Logan頭模作為實驗對象. 在0°～180°范圍均勻采集60個角度，探測器個數(shù)為260個. 實驗中，在濾波反投影(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

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

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

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

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

(23)

(24)

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

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

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

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

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

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

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

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

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

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

5 結 論

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