結(jié)合鄰域信息的TV正則化稀疏角度重建算法

齊澤瑤，王遠(yuǎn)軍

(上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海200093)

E-mail:yjusst@126.com

1 引言

計算機(jī)斷層成像(Computed Tomography，CT)是醫(yī)學(xué)影像診斷領(lǐng)域中廣泛應(yīng)用的一種現(xiàn)代成像技術(shù).然而在醫(yī)療診斷過程中過度的輻射劑量可能會影響人體的生理機(jī)能、增加致癌機(jī)率，尤其是對兒童影響更為嚴(yán)重［1］.因此，在臨床檢查中輻射劑量的控制非常重要.在臨床檢查過程中劑量的高低控制主要通過X射線的強(qiáng)度和掃描投影的角度數(shù)量來調(diào)節(jié).當(dāng)X射線的強(qiáng)度減小時，采集的投影數(shù)據(jù)會不可避免地引入不穩(wěn)定噪聲，將會導(dǎo)致圖像質(zhì)量下降.而減少投影角度的數(shù)量屬于不完全投影數(shù)據(jù)重建問題，是數(shù)學(xué)理論中的一個病態(tài)逆問題［2］.而壓縮感知(Compressed Sensing，CS)理論的發(fā)展為不完全投影重建提供了新思路:如果圖像是稀疏的，或可以找到合適的變換域使其稀疏，則可以通過挖掘其稀疏圖像的L1范數(shù)準(zhǔn)確地重建圖像［3］.

在欠采樣稀疏投影數(shù)據(jù)的條件下重建圖像時，由于投影信息的不足，采用解析算法，如濾波反投影(filtered back-projection，F(xiàn)BP)算法時會產(chǎn)生混疊偽影，這將極大地影響重建圖像的質(zhì)量;采用傳統(tǒng)的迭代重建算法，如代數(shù)重建算法(algebraic reconstruction technique，ART)和聯(lián)合代數(shù)重建算法(simultaneous algebraic reconstruction technique，SART)雖然不會受到混疊偽影的干擾可以重建出較清晰的圖像，但也會產(chǎn)生一些不符合臨床需要的偽影從而無法獲得理想的結(jié)果;而依據(jù)CS理論通過在目標(biāo)函數(shù)中加入懲罰項或正則項時可以進(jìn)一步改善重建圖像的質(zhì)量，其中以全變分(total variation，TV)作為CS稀疏變換的重建算法取得了很好的重建效果并得到了廣泛的研究［4-8］.

2006年，Sidky等［9］將TV模型應(yīng)用到了圖像重建領(lǐng)域，提出了TVM-POCS(total variation minimization-projection onto convex sets)算法，表明了以TV為先驗信息作為正則化可以在稀疏投影數(shù)據(jù)中重建出高質(zhì)量的圖像，但由于圖像的分段常數(shù)假設(shè)，在重建結(jié)果圖像的邊緣和細(xì)節(jié)區(qū)域有時會過于平滑.為了解決這個問題，一些研究學(xué)者考慮到各向同性邊緣屬性的局限性，更充分的利用部分相關(guān)先驗信息，以權(quán)重的形式引入到TV項來利用各向異性邊緣屬性提高重建算法的性能，并提出相關(guān)的算法通過在圖像的邊緣和平滑區(qū)域引入自適應(yīng)權(quán)重保護(hù)了圖像的邊緣和細(xì)節(jié)信息［10-14］.

Tian等［10］提出 EPTV(TV-based edge preserving)模型，通過局部梯度強(qiáng)度構(gòu)建自適應(yīng)的指數(shù)權(quán)重函數(shù)并引入到TV的每一項中，保護(hù)了更多的低對比度結(jié)構(gòu)和邊緣信息.然而，僅考慮了各向同性邊緣屬性;Liu等［11］提出AwTV(adaptiveweighted TV)模型，通過局部梯度強(qiáng)度構(gòu)建自適應(yīng)的指數(shù)權(quán)重函數(shù)引入到TV模型中的各個梯度的子方向上，考慮了各向異性邊緣屬性以圖像局部信息自適應(yīng)調(diào)整梯度強(qiáng)度，更好的保存了邊緣細(xì)節(jié)信息.然而，也僅僅考慮了圖像的局部信息;Zhang等［12］提出了 NLTV(Non-local TV)模型，通過利用圖像的非局部自相似性特質(zhì)構(gòu)建自適應(yīng)的權(quán)重函數(shù)引入到TV模型中，更好的保存了邊緣細(xì)節(jié).綜合考慮局部范圍和遠(yuǎn)程范圍的相關(guān)性，范圍的大小與重建時間成正比，當(dāng)局部范圍較大時耗時較長;Rong等［13］在AwTV模型的基礎(chǔ)上提出了EGAwTV(anisotropic edge-guided adaptive-weighted TV)模型，通過檢測邊緣信息和局部圖像強(qiáng)度梯度來構(gòu)建權(quán)重函數(shù)，引入更多的先驗信息進(jìn)一步保存了圖像的邊緣和細(xì)節(jié)信息.然而，僅僅考慮了圖像的一階信息;Zheng等［14］提出了TACTV(total absolute curvature TV)模型，將曲率的絕對值引入到TV模型中，利用圖像的二階信息能夠更好地描述特征復(fù)雜的圖像等等.

受到上面提及思想的啟發(fā)，本文在傳統(tǒng)的TV重建算法的基礎(chǔ)上引入一個權(quán)重來自適應(yīng)的調(diào)節(jié)橫向梯度和縱向梯度在TV中的占有比.其中的權(quán)重大小取決于每個像素鄰域區(qū)域內(nèi)的均值和均方差.通過像素鄰域信息的均值和均方差構(gòu)建而成的自適應(yīng)權(quán)重利用各向異性邊緣屬性來改進(jìn)傳統(tǒng)TV重建算法的性能.

2 本文算法

2.1 TV重建算法

CT圖像重建的數(shù)學(xué)模型通常采用線性方程組來描述，模型的離散形式表示為:

其中，p為投影向量，f為待重建圖像向量，系統(tǒng)矩陣A關(guān)聯(lián)了圖像向量f和投影向量p，其元素通常被設(shè)置為射線路徑與路徑上相關(guān)像素的相交長度［15］.CT圖像重建的問題可以表達(dá)為根據(jù)系統(tǒng)矩陣A和投影向量p來確定未知圖像向量f.2004年，Candes等提出的CS理論證明了稀疏信號可以從遠(yuǎn)小于奈奎斯特采樣要求的采樣數(shù)據(jù)中被重建.因此，在迭代重建算法求解時，可以應(yīng)用包括圖像稀疏性的一些先驗信息［16］.CS理論在CT重建中的應(yīng)用可以將重建問題表示為帶約束的優(yōu)化問題:

對于醫(yī)學(xué)圖像來說大部分是非稀疏的，而圖像中某些結(jié)構(gòu)的邊緣處經(jīng)常發(fā)生急劇的強(qiáng)度變化，可以通過離散梯度變換進(jìn)行稀疏表示［3］.因此，稀疏投影數(shù)據(jù)重建經(jīng)常采用TV(圖像梯度的L1范數(shù))作為目標(biāo)函數(shù)進(jìn)行圖像的優(yōu)化［9］: 表示圖像在(s，t)處的像素值，s和t分別表示橫向和縱向坐標(biāo).

2.2 本文算法

為了進(jìn)一步改善TV重建算法的性能，本文通過聯(lián)合圖像像素鄰域信息的均值和均方差利用各向異性邊緣屬性提出了一種結(jié)合鄰域信息的TV正則化稀疏角度重建算法.自適應(yīng)權(quán)重函數(shù)構(gòu)建原則如下［17］:遍歷原圖像，以每個像素點(diǎn)為中心，計算其鄰域內(nèi)像素的均值和均方差，然后計算出權(quán)值系數(shù) ws，t:

其中，

自適應(yīng)權(quán)重 ws，s－1，t，t和 ws，s，t，t－1分別以橫縱方向上的梯度圖像為基礎(chǔ)，采用公式(4)構(gòu)建而成.

TV重建算法將圖像作為一個整體考慮，有效的利用了圖像的稀疏性和邊緣方向先驗信息，能夠在稀疏投影采樣的條件下重建出質(zhì)量較高的圖像，然而由于各向同性邊緣屬性的局限性不能充分的調(diào)整局部信息，重建圖像往往會出現(xiàn)邊緣過于平滑的現(xiàn)象.而本文提出的重建算法采用權(quán)重函數(shù)自適應(yīng)的調(diào)節(jié)橫向梯度和縱向梯度在TV中的占有比利用各向異性邊緣屬性打破了局限性.

圖像中每個像素在相同鄰域大小的區(qū)域內(nèi)分布都是不同的，這里我們給定一個系數(shù)w來表示它們之間的不同.而均值和均方差可以較好的體現(xiàn)各個像素鄰域區(qū)域內(nèi)像素的分布以及變化程度.均值表示鄰域區(qū)域內(nèi)像素整體的分布水平，均方差表示大部分像素值和其均值的偏離程度反映了鄰域區(qū)域內(nèi)像素波動大小的程度.當(dāng)鄰域區(qū)域內(nèi)像素值大小越接近時說明此區(qū)域越趨于平滑，均方差越小，反之均方差越大.這兩者共同反映了圖像像素鄰域區(qū)域的系數(shù)w的大小，從而在不同像素鄰域區(qū)域內(nèi)進(jìn)行不同尺度的正則化實(shí)現(xiàn)了局部信息的調(diào)整.

2.3 本文算法實(shí)現(xiàn)及步驟

受到傳統(tǒng)TV重建算法實(shí)現(xiàn)的啟發(fā)，本文算法實(shí)現(xiàn)也采用最速下降法求解此優(yōu)化問題.當(dāng)圖像質(zhì)量遠(yuǎn)遠(yuǎn)達(dá)不到理想的效果時，圖像像素鄰域區(qū)域的均值和均方差會受到很大的影響.本文算法的實(shí)現(xiàn)過程分為兩部分:在執(zhí)行梯度下降算法選擇下降方向優(yōu)化圖像時，首先使用傳統(tǒng)的TV作為下降方向，當(dāng)?shù)欢ǖ拇螖?shù)后再選擇帶權(quán)重的TV作為下降方向.本文算法實(shí)現(xiàn)的主要步驟如下:

1)初始化重建參數(shù):

其中，n≤N是算法當(dāng)前迭代的次數(shù)，N是總的迭代次數(shù).

2)利用ART算法重建過渡圖像:

其中，λ 是松弛因子，m=1，…，Ndata，Ndata是所有投影角度總的射線條數(shù).

3)對重建圖像進(jìn)行非負(fù)性約束:

4)選擇像素鄰域區(qū)域大小，計算權(quán)重:

5)利用梯度下降算法優(yōu)化圖像，其中l(wèi)=1，2，…，M:

5-3)更新重建圖像f:

其中，M為梯度下降算法的迭代次數(shù)，k是算法迭代TV的總次數(shù)，ξ是為了防止分母為零而設(shè)置的一個很小的正數(shù)，這里設(shè)置為ξ=10－8.

6)循環(huán)2)-5)，直到滿足收斂標(biāo)準(zhǔn)或得到最大迭代次數(shù)為止，本文中我們設(shè)置達(dá)到迭代次數(shù)N為停止迭代標(biāo)準(zhǔn).

3 數(shù)值仿真

為了研究本文算法在欠采樣稀疏投影采樣的條件下重建的可行性，我們在此章節(jié)對仿真數(shù)據(jù)Shepp-Logan體模［18］和真實(shí)的核桃投影數(shù)據(jù)［19］進(jìn)行重建實(shí)驗.所有實(shí)驗測試使用的電腦配置為 Intel(R)Core(TM)i7-7700 CPU，主頻為3.60GHz，運(yùn)行內(nèi)存為 8GB，操作系統(tǒng)為 Windows 8，測試平臺為MATLAB 2016b.

3.1 Shepp-Logan體模仿真實(shí)驗

本章節(jié)實(shí)驗以Shepp-Logan模型模擬生成基于圓形軌道的扇形束投影數(shù)據(jù)，仿真模型的大小設(shè)為256*256.在模擬環(huán)境中，使用等角虛擬探測器，扇角設(shè)為π/4，射線源與旋轉(zhuǎn)中心的距離設(shè)為512mm，探測器通道數(shù)設(shè)置為256個.在2π范圍內(nèi)對體模進(jìn)行均勻采樣，分別獲取24、36、48個投影數(shù)據(jù).重建過程中的重建參數(shù)設(shè)置為:總的迭代次數(shù)N=200，ART迭代算法的松弛因子λ=0.25，像素鄰域區(qū)域的大小p*q=3*3，TV 的迭代次數(shù)為 k=50，α =0.2，M=20.圖 1 為TV算法和本文算法重建得到的結(jié)果圖像.

圖1 TV算法(第1行)和本文算法(第2行)從24(第1列)、36(第2列)、48(第3列)個無噪聲投影數(shù)據(jù)重建的結(jié)果圖像;第3行和第4行分別是TV算法和本文算法結(jié)果圖像的感興趣區(qū)局部放大顯示Fig.1 Reconstructed results of the TV algorithm(1st line)and our algorithm(2nd line)from 24(1st column)，36(2nd column)and 48(3rd column)noise-free projection data.The 3rd line and 4th line are respectively the corresponding zoomed in part display of ROI of the result image of the TV algorithm and our algorithm

從圖1中的第1行和第2行我們可以看出兩種迭代算法在不同的投影數(shù)據(jù)下重建圖像都基本無偽影，但在投影數(shù)據(jù)較少時的重建圖像在邊緣的細(xì)節(jié)區(qū)域顯示較模糊，而投影數(shù)據(jù)較多時的重建圖像在邊緣的細(xì)節(jié)區(qū)域顯示較清晰，在視覺效果上與原始的Shepp-Logan體模基本無差異.為了能夠更好的比較兩種算法，部分感興趣區(qū)域已經(jīng)被白色虛線框標(biāo)出并且放大顯示在圖1的第3行和第4行中，從中可以清晰的發(fā)現(xiàn)，本文算法無論是在投影數(shù)據(jù)較少時還是在投影數(shù)據(jù)較多時，重建的結(jié)果圖像邊緣都比TV算法重建的結(jié)果圖像清晰.因此，本文算法可以重建出質(zhì)量高于TV算法的圖像.為了進(jìn)一步驗證本文算法的優(yōu)越性能，利用均方根誤差(root mean square error，RMSE)表明算法的重建精度.RMSE定義如下:

其中，fj和gj分別表示原圖像向量和重建圖像向量，RMSE通常用來評估重建圖像的精確度，它的取值范圍為0-1，其值越接近0表明重建精度越高.

圖2給出了兩種重建算法分別在24、36和48個無噪聲投影數(shù)據(jù)下重建結(jié)果關(guān)于RMSE隨迭代次數(shù)增加的變化曲線圖.本文算法前k=50次迭代選擇TV下降方向進(jìn)行迭代，因此前面RMSE下降趨勢相同.但在圖2中可以清晰看出，本文算法比TV重建算法在相同的迭代次數(shù)時重建結(jié)果誤差較小.綜合前面的分析，可以得出本文算法具有比TV重建算法更加卓越的重建性能.

圖2 TV算法和本文算法在24(第1列)、36(第2列)、48(第3列)個無噪聲投影數(shù)據(jù)下重建結(jié)果關(guān)于RMSE隨迭代次數(shù)增加的變化曲線Fig.2 Under the condition of 24(1st column)，36(2nd columns)and 48(3rd columns)noise-free projection data，the RMSE curves of the reconstructed image of TV algorithm and our algorithm changes with the increase of iterations

3.2 算法抗噪聲性能的驗證

為了驗證算法的抗噪聲性能，在獲得的Shepp-Logan模型投影數(shù)據(jù)中加入一個隨機(jī)噪聲，然后分別降采樣為24、36和48個含有噪聲的投影數(shù)據(jù)，采用TV重建算法和本文算法進(jìn)行重建實(shí)驗.重建過程中的重建參數(shù)設(shè)置與2.1節(jié)一致.兩種重建算法的結(jié)果圖像和部分感興趣區(qū)域局部放大部分顯示在圖3中.

圖3 TV算法(第1行)和本文算法(第2行)從24(第1列)、36(第2列)、48(第3列)個含噪聲投影數(shù)據(jù)重建的結(jié)果圖像;第3行和第4行分別是TV算法和本文算法結(jié)果圖像的感興趣區(qū)局部放大顯示Fig.3 Reconstructed results of the TV algorithm(1st line)and our algorithm(2nd line)from 24(1st column)，36(2nd column)and 48(3rd column)noise projection data.The 3rd line and 4th line are respectively the corresponding zoomed in part display of ROI of the result image of the TV algorithm and our algorithm

從圖3中我們可以看出重建圖像基本沒有偽影，在相同迭代次數(shù)下重建圖像的邊緣區(qū)域隨著投影數(shù)據(jù)的增加而逐漸的清晰.對比圖1和圖3可以從感興趣區(qū)域看出，在邊緣處不如無噪聲投影數(shù)據(jù)重建結(jié)果清晰.從圖3中也可以清晰的看出相同數(shù)量的投影數(shù)據(jù)重建結(jié)果中，本文算法重建的結(jié)果圖像比TV重建算法重建的結(jié)果圖像邊緣處較清晰.圖4給出了兩種重建算法分別在24、36、48個含噪聲投影數(shù)據(jù)的條件下重建的結(jié)果關(guān)于RMSE隨迭代次數(shù)增加的變化曲線.可以看出本文算法比TV重建算法更快的收斂，且重建結(jié)果誤差較小.這與在無噪聲的投影數(shù)據(jù)重建實(shí)驗結(jié)果一致.因此，我們可以得出本文算法比TV算法具有更好的重建性能和更快的收斂速度.

圖4 TV算法和本文算法在24(第1列)、36(第2列)和48(第3列)個含噪聲投影數(shù)據(jù)下重建結(jié)果關(guān)于RMSE隨迭代次數(shù)增加的變化曲線Fig.4 Under the condition of 24(1st column)，36(2nd columns)and 48(3rd columns)noise projection data，the RMSE curves of the reconstructed image of TV algorithm and our algorithm changes with the increase of iteration

3.3 真實(shí)的投影數(shù)據(jù)重建

為了更進(jìn)一步驗證算法性能，本文對獲得的真實(shí)CT投影數(shù)據(jù)進(jìn)行稀疏角度重建實(shí)驗.被檢測的物體為一個核桃模型，投影數(shù)據(jù)是用德國菲尼克斯(Phoenix|X-Ray)X射線檢測系統(tǒng)與服務(wù)有限公司提供專有的μCT成像系統(tǒng)掃描獲得.其中，射線源到平板探測器的距離為300mm，射線源到掃描物體的距離為110mm，平板探測器的寬度為114.8mm.在管電壓為80kV，管電流為200μA時的狀態(tài)下以圓軌跡錐束掃描進(jìn)行投影數(shù)據(jù)的獲取.然后降采樣分別獲得30、60、120個角度的投影數(shù)據(jù).

圖5 TV算法(第1行)和本文算法(第2行)從30(第1列)、60(第2列)、120(第3列)個投影數(shù)據(jù)重建的結(jié)果圖像Fig.5 Reconstructed results of the TV algorithm(1st line)and our algorithm(2nd line)from 30(1st column)，60(2nd column)and 120(3rd column)projection data

本文對z軸上的中心切片進(jìn)行重建，TV算法和本文算法重建參數(shù)設(shè)置與仿真實(shí)驗一致，重建結(jié)果和部分感興趣區(qū)局部放大如圖5和圖6所示.從圖5中可以看出重建的結(jié)果圖像隨著投影數(shù)據(jù)的增加，雖然邊緣有模糊的現(xiàn)象但在相同迭代次數(shù)下圖像質(zhì)量越來越好.結(jié)合圖6顯示的部分感興趣區(qū)的局部放大圖可知在視覺效果上本文算法的重建結(jié)果較優(yōu).

圖6 TV算法(第1行)和本文算法(第2行)從30(第1列)、60(第2列)和120(第3列)個投影數(shù)據(jù)重建結(jié)果圖像感興趣區(qū)Fig.6 Region of interest of reconstructed results of the TV algorithm(1st line)and our algorithm(2nd line)from 30(1st column)，60(2nd column)and 120(3rd column)projection data

4 結(jié)論

本文在稀疏角度投影數(shù)據(jù)的條件下，針對傳統(tǒng)的TV算法的局限性提出了一種結(jié)合鄰域信息的TV正則化稀疏角度重建算法.與TV重建算法相比，本文算法通過像素鄰域區(qū)域的均值和均方差引入了一個權(quán)重，考慮了圖像的各向異性邊緣屬性，使圖像的局部信息可以進(jìn)行自適應(yīng)的調(diào)整，并運(yùn)用最速下降法進(jìn)行求解，更好的提高了重建圖像的質(zhì)量.并且以Shepp-Logan仿真模型和核桃的真實(shí)投影數(shù)據(jù)進(jìn)行了重建實(shí)驗，實(shí)驗結(jié)果也表明了本文算法的優(yōu)越性能，可以在抑制偽影的同時更好地保留邊緣結(jié)構(gòu)信息，重建更精確的圖像，驗證了本文提出的算法的有效性.

本文只對仿真實(shí)驗和核桃中心層面的圖像進(jìn)行了重建，實(shí)驗數(shù)據(jù)有限，有待進(jìn)一步深入研究.另外，本文重建算法采用梯度下降法求解，耗時長，可以進(jìn)一步改進(jìn)現(xiàn)有技術(shù)，如基于GPU的加速技術(shù)和新的求解算法.此外，我們將嘗試探索有限角度投影重建問題，以減少輻射劑量.