莫禮東，徐伯慶

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

?

一種改進(jìn)的TV圖像重建算法

莫禮東，徐伯慶

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

針對(duì)當(dāng)前稀疏角度下有限角圖像重建過程中，邊界部分出現(xiàn)偽影，降低了圖像重建質(zhì)量的缺陷。文中提出了一種新的ART+TV算法，該方法是在原始TV算法的基礎(chǔ)上進(jìn)行改進(jìn)。原始TV梯度下降算法求解目標(biāo)函數(shù)最小值時(shí)，使用固定函數(shù)作為目標(biāo)函數(shù)，文中對(duì)其進(jìn)行更改，采用帶參數(shù)的目標(biāo)函數(shù)，并對(duì)TV重建后的結(jié)果進(jìn)行自適應(yīng)步長(zhǎng)修正，加速圖像收斂。與傳統(tǒng)的ART+TV算法相比，文中算法在不改變重建速度的基礎(chǔ)上，且在少量迭代次數(shù)下，能重建出質(zhì)量更高的圖像，抑制圖像偽影。

ART+TV算法；迭代重建；稀疏角度CT

MO Lidong, XU Boqing

(School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology,Shanghai 200093, China)

X射線CT成像已在醫(yī)學(xué)和工業(yè)中有廣泛應(yīng)用[1]。但由于X射線的輻射劑量與X射線曝光次數(shù)成正比，而X射線的輻射對(duì)人體的身體健康有較大傷害。因此，減少X射線的輻射劑量成為了醫(yī)學(xué)圖像重建的研究方向。通過控制CT系統(tǒng)。對(duì)X射線投影角度進(jìn)行稀疏采樣和有限角度下的CT圖像重建算法，是一種有效減少病人輻射劑量的技術(shù)手段，并被廣泛應(yīng)用。本文主要研究有限角度下的圖像重建問題。

傳統(tǒng)的圖像重建算法主要包括濾波反投影(FBP)算法和代數(shù)迭代(ART)算法。其中FBP算法是從連續(xù)的成像模型推導(dǎo)出來的，并且間接認(rèn)為在掃描角度范圍內(nèi)進(jìn)行高密度的投影采樣[2]。因此，該方法難以實(shí)現(xiàn)有限角度重建算法。ART算法是代數(shù)迭代重建，相比FBP算法，此算法更加實(shí)用于稀疏角度圖像重建。文獻(xiàn)[3～4]提出一種基于最小化和投影到凸集的迭代重建算法，且應(yīng)用到稀疏角度CT重建，并取得了較好的重建結(jié)果。傳統(tǒng)的ART+TV算法在掃描角度比較少的情況下，重建結(jié)果出現(xiàn)變形和偽影現(xiàn)象。為客服這一缺點(diǎn)，不少算法在原算法的基礎(chǔ)上進(jìn)行改進(jìn)。文獻(xiàn)[5～6]提出了一種基于先驗(yàn)圖像約束的壓縮感知重建算法。其主要思想是通過完整數(shù)據(jù)重建出來的圖像作為先驗(yàn)條件，然后獲取多個(gè)角度下的投影數(shù)據(jù)，將TV正則化思想用到待重建圖像的基礎(chǔ)上，還用到待重建圖像與先驗(yàn)圖像的差值圖像。能夠用少量的二維投影數(shù)據(jù)進(jìn)行高質(zhì)量的圖像重建。但此算法對(duì)先驗(yàn)圖像的信息具有一定的條件，給算法帶來了局限性。本文在傳統(tǒng)ART+TV算法的基礎(chǔ)上，以ART算法為基礎(chǔ)，加上全變差約束條件，對(duì)梯度圖像的計(jì)算進(jìn)行修改，采用可變的梯度計(jì)算函數(shù)，在保證重建圖像速度的基礎(chǔ)上，進(jìn)行質(zhì)量的優(yōu)化。通過大量實(shí)驗(yàn)結(jié)果驗(yàn)證本文算法在原始算法的基礎(chǔ)上確實(shí)有所改進(jìn)。

1 ART+TV算法原理

1.1 ART算法原理

ART迭代算法最初是于1937年由Kaczmarz提出的[7-8]，其旨在求解類似式(1)的方程組

(1)

在圖像重建中，也是一個(gè)大稀疏矩陣的求解。因此，方程(1)的求解與圖像重建求解方法類似。

假設(shè)一條強(qiáng)度為I0的單射線穿過一個(gè)內(nèi)部分布未知的物體后，檢測(cè)器上檢測(cè)到的射線的強(qiáng)度為I，射線在物體中的衰減過程滿足比爾定律[9]即

(2)

其中，f(x,y)表示物體在x-y平面內(nèi)對(duì)單射線L的線性衰減系數(shù)。由式(2)可看出,f的值越大，衰減程度就越大。令穿過物體后射線投影數(shù)據(jù)為p，則有

(3)

因此p對(duì)應(yīng)f(x,y)沿射線L方向的線積分。ART算法首先將待重建圖像分解成n×n的小格，每一個(gè)小格對(duì)應(yīng)一個(gè)像素值，如圖1所示。

圖1 ART重建示意圖

在圖1中,fi所表示的是第i個(gè)像素值，i=1,2,...,N,N=n2表示重建圖像包含的總的像素?cái)?shù)，因此f=(f1,f2,...,fN)表示一個(gè)N維空間的向量。在此我們定義掃描射線的寬度與像素網(wǎng)格的寬度一致。若在掃描時(shí)有P個(gè)方向，每個(gè)方向有R條射線，則總的射線條數(shù) 。圖像像素在第j條射線上的積分記做Pj,j=1,2,...,M,所以圖像向量f和投影向量p的關(guān)系可表示為

(4)

式(4)中，wij為加權(quán)因子，也即系統(tǒng)矩陣，反映了第j個(gè)像素對(duì)第i條射線投影值的貢獻(xiàn)。每一條掃描射線的投影過程即為上述的一個(gè)方程，解得線性方程組的解也就得到了整個(gè)待重建截面的像素密度分布。

而由系數(shù)wij構(gòu)成的矩陣是一個(gè)大型的稀疏矩陣，用常規(guī)的方法計(jì)算不僅計(jì)算過程復(fù)雜，且計(jì)算速度緩慢。在實(shí)際中，均用迭代方法進(jìn)行計(jì)算，也即ART算法，其迭代公式為

(5)

式(5)中，k為迭代次數(shù)，wij為投影系數(shù)。由于噪聲的影響，式(5)求解的結(jié)果就會(huì)存在較大的誤差，于是引入松弛因子，來減少重建過程中的誤差，其改進(jìn)公式為

(6)

1.2 TV原理

在醫(yī)療圖像重建中，由于X射線對(duì)人體的傷害以及人體器官組織的相似性。因此，稀疏角度圖像重建被廣泛運(yùn)用?；趬嚎s感知的原理，可由少量的投影數(shù)據(jù)重建出原始圖像，其被稱為ART+TV算法。其重建圖模型可表示為

min‖f‖TV,s.t.Wf=p

(7)

式(7)中,W是投影系數(shù)矩陣;p是投影值矩陣;f是灰度化的像素值。式(7)與式(1)是對(duì)應(yīng)的，只是在式(1)的基礎(chǔ)上加入了約束條件，其也是ART+TV與ART的主要區(qū)別。由式(7)可知，求解min‖f‖TV是關(guān)鍵。而并不能直接求出min‖f‖TV，一般均是通過用1范數(shù)L1來求解。

式(8)中▽fx,y=(Dxf,Dyf)，Dxf和Dyf是離散微分算子，對(duì)于圖像上的點(diǎn)(x,y)，其梯度計(jì)算跟f(x,y)和f(x,y)附近的像素值有關(guān)

Dxfx,y=fx,y-fx-1,y;Dyfx,y=fx,y-fx,y-1

(8)

L1范數(shù)就是圖像中每一個(gè)像素值的絕對(duì)值之和，其求導(dǎo)運(yùn)算不方便。因此，大多采用L2范數(shù)來近似L1范數(shù)即

(9)

因此，對(duì)某個(gè)像素值f(x,y)求導(dǎo)，得

式中，ε=10-8其主要是防止‖f‖TV求導(dǎo)的時(shí)候出現(xiàn)分母為零的情況。通過上式的求解步驟，就容易求解方程(7)。

傳統(tǒng)的ART+TV算法步驟如下[10]

Givenpi,Nview,Wi(i=1,2,3,…Nview)andNTVInitialization:f0=0.

fork=1,2,3,…do

ART Updating:

f0=fk

fori=1,2,3,…Nviewdo

end

TV Minimization:

forn=1,2,3,…,NTV-1 do

d=‖fn-f0‖2,

fn+1=fn+a·d·v

end

其中，Nview是ART迭代次數(shù);NTV是TV的迭代次數(shù);pi是投影值;Wi是投影系統(tǒng)矩陣，本實(shí)驗(yàn)中NTV取20，Nview取10，K取20，a取0.2。

2 本文算法

本文是在原始的ART+TV的算法上進(jìn)行改進(jìn)。原始的ART+TV算法，當(dāng)用L1范數(shù)計(jì)算其最小值時(shí)，對(duì)于邊界的方向，無法檢測(cè)出來，由式(8)可知，L1范數(shù)求解最小值是求和的形式求解的，其中包含了方向不確定的邊界的貢獻(xiàn)，這樣對(duì)重建圖像會(huì)產(chǎn)生偽影。因此，本文提出了一種新的TV算法，對(duì)式(3)進(jìn)行改進(jìn)，采用式(10)進(jìn)行替換

(10)

因此，此時(shí)對(duì)‖f‖TV求導(dǎo)公式變?yōu)?/p>

并對(duì)傳統(tǒng)梯度下降的固定步長(zhǎng)d進(jìn)行修正，加速圖像的收斂，其改動(dòng)原理為

d2=‖fNview-fn+1‖

(11)

若d2>0.85×d,則d=0.85×d2。通過修改后的實(shí)驗(yàn)步驟為：(1)初始化圖像及相關(guān)參數(shù)；(2)根據(jù)已經(jīng)投影數(shù)據(jù)進(jìn)行ART迭代，并進(jìn)行非負(fù)修正，得到圖像1；(3)利用本文提出的帶參數(shù)求解的方法進(jìn)行TV梯度下降算法求解目標(biāo)函數(shù)最小值，得到重建圖像2；(4)并對(duì)傳統(tǒng)的TV算法進(jìn)行動(dòng)態(tài)步長(zhǎng)d2修正，如式(11)所示，加速圖像收斂，得圖像3；(5)判斷是否滿足迭代終止條件，則退出循環(huán)；否則返回步驟(2)，將圖像3作為步驟(2)中ART迭代的初始圖像，循環(huán)步驟(2)～(5)。

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

本實(shí)驗(yàn)采用分辨率為128×128的經(jīng)典Shepp-Logan模型進(jìn)行實(shí)驗(yàn)，分別在0～90；0～120；0～150范圍，間隔1度進(jìn)行仿真重建，對(duì)本文算法，傳統(tǒng)ART+TV算法，ART算法進(jìn)行對(duì)比，以說明改進(jìn)后的算法的優(yōu)勢(shì)，其實(shí)驗(yàn)結(jié)果如圖2所示。

圖2 3組實(shí)驗(yàn)結(jié)果圖

如圖2所示，第1組圖像實(shí)驗(yàn)是直接進(jìn)行ART代數(shù)迭代算法的實(shí)驗(yàn)結(jié)果圖；第2組圖像實(shí)驗(yàn)是ART+TV即傳統(tǒng)的TV算法實(shí)驗(yàn)；第3組是本文算法的實(shí)驗(yàn)結(jié)果圖。

在實(shí)驗(yàn)過程中，為了客觀地評(píng)價(jià)出重建圖像的質(zhì)量，本文采用峰值信噪比(PSNR)歸一化均方距離判距d來定量評(píng)價(jià)重建圖像與原始圖像的誤差。表1顯示了3組實(shí)驗(yàn)的參數(shù)指標(biāo)對(duì)比，其中參數(shù)PSNR值越大說明圖像噪聲越小，圖像質(zhì)量越好；歸一化均方距離判距d的值越小代表重建誤差越小。

表1 90個(gè)投影角度下重建指標(biāo)對(duì)比

表2 120個(gè)投影角度下重建指標(biāo)對(duì)比

表3 150個(gè)投影角度下重建指標(biāo)對(duì)比

從表1～表3的實(shí)驗(yàn)結(jié)果可看出，本文算法在傳統(tǒng)算法的基礎(chǔ)上有所改進(jìn)。

圖3是在0°～90°范圍內(nèi)不同的迭代次數(shù)下的PSNR值對(duì)比圖。

圖3 不同算法的PSNR值的對(duì)比圖

如圖3可示，本文算法的PSNR值比其他兩種算法的值要大。所以，本文算法對(duì)有限角度下圖像重建算法的重建圖像質(zhì)量有所提高。

4 結(jié)束語

在稀疏角度圖像重建中，本文針對(duì)有限角度圖像

重建進(jìn)行了研究。在有限角度重建過程中，受重建角度的限制，重建角度少的情況下邊界出現(xiàn)偽影，為了避免這種情況，基于ART+TV算法提出了一種算法，其對(duì)梯度圖像的求解過程進(jìn)行修正，并且在梯度下降算法中加入動(dòng)態(tài)修正因子，加速收斂。通過實(shí)驗(yàn)驗(yàn)證，本文算法在不改變重建的速度的基礎(chǔ)上，比ART+TV算法在圖像重建的質(zhì)量上有所改進(jìn)。

[1] 余維.不完備投影數(shù)據(jù)的CT重建算法研究[D].重慶:重慶大學(xué),2014.

[2] Pan X,Sidky E Y,Vannier M.Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction[J].Inverse Problems,2009,25(12):1541-1548.

[3] Gopi V P,Fayiz T K,Palanisamy P. Regularization based CT image reconstruction using algebraic techniques[C].Washington,USA:International Conference on Electronics and Communication Systems, IEEE,2014.

[4] Sidky E Y,Kao C M,Pan X.Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J].Journal of X-ray science and technology,2009,14(2):119-139.

[5] Chen G J,Leng S.Prior image constrained compressed sensing (PICCS):a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets[J].Medical Physics,2008,35(2):660-663.

[6] 郭靜鈺,齊宏亮,袁媛,等.先驗(yàn)圖像約束的有限角度CT圖像重建算法[J].核電子學(xué)與探測(cè)技術(shù),2014(12):1421-1424.

[7] Kaczmarz S.Angenaherte auflosung von systemen lincarer glerichungen[J].Bull Acad P0l Sciencever Letter A,1937(6-8A):355-357.

[8] Tanabe K.Projection method for solving a singular system of linear equations and its applications[J].Numer Math,1971,17(3):203-214.

[9] G T赫爾曼.由投影重建圖像-CT的理論基礎(chǔ)[M].北京:科學(xué)出版社,1985.

[10] Jin X,Li L,Chen Z,et al.Anisotropic total variation for limited-angle CT reconstruction[C].Tianjin:IEEE Nuclear Science Symposium Conference Record, Nuclear Science Symposium,2010.

An Improved TV Image Reconstruction Algorithm

For the current limited angle image reconstruction process sparse angle, the boundary part artifacts appear to reduce the image quality of the reconstructed defects. This paper proposes a new algorithm of ART + TV, which is the basis of the original TV algorithms, to improve it. In solving the minimum objective function by the original TV gradient descent algorithm, the fixed function is used as the objective function. While in this paper, the text to change it, the objective function with parameters is employed, and the results of TV reconstruction are adaptive step corrected to accelerate the graphics convergence. With compare the proposed algorithm offers better quality of reconstructed images with a few iterations than the conventional ART + TV algorithms at the same reconstruction speed and suppresses image artifacts

ART+TV algorithm; iterative reconstruction; sparse angle CT

2016- 01- 14

莫禮東(1991-)，男，碩士研究生。研究方向：信號(hào)與信息處理。徐伯慶(1958-)，男，博士，副教授。研究方向：通信及圖像處理。

10.16180/j.cnki.issn1007-7820.2016.10.014

TP391.41

A

1007-7820(2016)10-047-04