各向異性邊緣檢測導(dǎo)引的CT稀疏角度重建

戎軍艷，劉文磊，高 鵬，廖琪梅，盧虹冰

各向異性邊緣檢測導(dǎo)引的CT稀疏角度重建

戎軍艷，劉文磊，高 鵬，廖琪梅，盧虹冰

目的：提出一種CT稀疏角度重建方法，用以降低CT的輻射劑量。方法：基于圖像各向異性的邊緣檢測信息，在圖像重建過程中更新圖像總變差各加權(quán)分量的加權(quán)權(quán)重，然后采用圖像總變差最小化方法重建圖像。結(jié)果：計(jì)算機(jī)模擬結(jié)果表明，與傳統(tǒng)的總變差最小化方法相比，該方法可以更好地重建圖像邊緣，且邊緣檢測的結(jié)果會(huì)影響重建效果。結(jié)論：采用該方法可以增強(qiáng)圖像邊緣，提高重建圖像的分辨率，改善重建圖像的質(zhì)量。

各向異性；邊緣檢測；總變差最小化；稀疏角度

0 引言

醫(yī)學(xué)研究表明，X線對(duì)人體有損傷，電離輻射可從分子和細(xì)胞水平對(duì)生物體造成損傷，照射X光可誘發(fā)癌癥、白血病或其他遺傳性疾病，因此，降低X線的輻射劑量對(duì)CT的臨床應(yīng)用具有重要的意義[1]。降低輻射劑量常用的方法有基于硬件改進(jìn)的方法、基于條件優(yōu)化的方法和基于算法改進(jìn)的方法[2]。在不改變硬件配置的條件下，減少投影角度是一種直接而有效的降低輻射劑量的方法。然而，減少投影角度會(huì)導(dǎo)致數(shù)據(jù)量小，圖像重建不能滿足Tuy-Smith數(shù)據(jù)完備性條件[3-4]。傳統(tǒng)的濾波反投影（filtered back-projection，F(xiàn)BP）類重建方法導(dǎo)致重建圖像分辨率下降、噪聲大、條形偽影嚴(yán)重[5]，不再適用。近幾年，利用信號(hào)的稀疏化將CT稀疏角度重建轉(zhuǎn)化為求解約束或無約束的目標(biāo)函數(shù)的最優(yōu)化問題的算法取得了顯著成果[6-8]。Sidky等[6-7]提出的凸集投影（projection onto convex sets）-總變差（total variation，TV）最小化（簡稱TV-POCS）方法采用梯度變換計(jì)算圖像的總變差，利用梯度下降法求解使得總變差最小化的重建圖像。TV-POCS方法用于離散扇束和錐束CT重建，結(jié)果顯示其相對(duì)于常用的最大似然-期望最大化方法（maximum likelihood expectation maximization，MLEM）與代數(shù)迭代法（algebraic reconstruction technique，ART）可得到更滿意的重建圖像。但是，該方法在計(jì)算圖像總變差時(shí)采用圖像像素梯度向量的模相加得到，未考慮圖像梯度的方向性，因此，在圖像求解過程中不可避免地帶來圖像邊緣的模糊和分辨率的下降。本文基于此方法，提出圖像各向異性邊緣檢測導(dǎo)引（anisotropic edge guided）的各向異性加權(quán)的TV最小化重建方法（簡稱Aniso-EGTV），采用各向異性加權(quán)的圖像總變差計(jì)算方法，利用圖像邊緣檢測信息更新權(quán)重，在一定程度上可提高圖像的分辨率。

1 算法原理

1.1 CT成像模型和TV-POCS重建方法

在迭代類重建算法中，通常將CT成像模型離散化。對(duì)于離散CT成像系統(tǒng)，其投影可用式（1）表示：

其中，p→為投影向量；M為投影矩陣，描述每個(gè)圖像像素點(diǎn)對(duì)各個(gè)探測器單元的貢獻(xiàn)權(quán)重；為待重建圖像，圖像上每個(gè)像素點(diǎn)用ui，j（i，j∈[0，L-1]）表示，圖像大小為L×L。

TV-POCS方法可表示為式（2）所描述的最優(yōu)化問題，即求滿足投影方程且使得圖像總變差最小化的u→。圖像的總變差用式（3）表示?？梢钥闯?，該方法采用圖像像素梯度向量的模相加計(jì)算圖像的總變差，而梯度向量的模為圖像像素沿x和y 2個(gè)方向的一階導(dǎo)數(shù)（即梯度的分量）的平方和相加后的平方根，其在x和y方向的一階導(dǎo)數(shù)的權(quán)重因子是相同的（即各向同性），均為1。未考慮圖像梯度的方向性，因此，在圖像求解過程中不可避免地帶來圖像邊緣的模糊和分辨率的下降。

1.2 Aniso-EGTV重建方法

將圖像總變差采用式（4）中各向異性加權(quán)的方式表達(dá)，gα代表中各個(gè)分量的權(quán)重因子，Dα表示對(duì)圖像取變分操作，表示像素（i，j）與其鄰域（k，l）變分的權(quán)重。計(jì)算中分量的個(gè)數(shù)由像素點(diǎn)的鄰域α決定，具體用N（l，j）表示像素（l，j）的鄰域，在此我們采用四鄰域計(jì)算。由此，類似于TVPOCS方法，各向異性加權(quán)的TV最小化方法用式（5）表示。

1.2.1 權(quán)重的更新

在四鄰域中，權(quán)重初始值設(shè)為π/4，我們將權(quán)重因子g構(gòu)造為L×L×4的三維矩陣，計(jì)算u→aniso-TV時(shí)每個(gè)像素點(diǎn)有4個(gè)權(quán)重因子與其對(duì)應(yīng)。權(quán)重隨著圖像迭代求解過程不斷更新，更新方法為：對(duì)選擇的圖像進(jìn)行邊緣檢測，當(dāng)圖像像素位于邊緣時(shí)將權(quán)重設(shè)為零，否則保持不變。

1.2.2 邊緣檢測方法

為了滿足權(quán)重因子的更新規(guī)則，我們改進(jìn)了傳統(tǒng)的Canny檢測方法[9]，提出各向異性的Canny邊緣檢測方法。具體過程如下：

（1）對(duì)圖像進(jìn)行高斯濾波。

（2）對(duì)濾波后的圖像分別沿x、y方向求導(dǎo)，得到對(duì)應(yīng)的梯度分量圖像Dx與Dy。

（3）沿梯度方向?qū)V波后的圖像進(jìn)行梯度模的非極大值抑制，確定梯度模的局部最大值點(diǎn)位置。

（4）對(duì)步驟（3）找到的梯度模的局部最大值點(diǎn)的梯度模進(jìn)行高、低雙閾值約束，從而找到大于等于高閾值的點(diǎn)域Ω和位于高、低閾值之間的點(diǎn)域Ψ。

（5）根據(jù)梯度分量圖像Dx與Dy將圖像像素點(diǎn)（i，j）按照其梯度矢量所屬區(qū)域分為x向邊緣點(diǎn)以及y向邊緣點(diǎn)。x向邊緣指邊緣切線方向偏向x軸，y向邊緣指邊緣切線方向偏向y軸。

（6）分別針對(duì)x向邊緣和y向邊緣，對(duì)Ω中的點(diǎn)以Ψ中的點(diǎn)作為連通鄰域補(bǔ)充進(jìn)行形態(tài)學(xué)連接和邊緣細(xì)化。

由此可得到沿x方向和y方向的2幅邊緣圖像，可利用這2幅邊緣圖像進(jìn)行上節(jié)中權(quán)重的更新。

1.2.3 Aniso-EGTV算法流程

參考TV-POCS算法，總結(jié)Aniso-EGTV流程如圖1所示，具體如下：

（1）設(shè)置初始迭代次數(shù)k=0，初始化重建圖像u（0），圖像大小為L×L；同時(shí)，初始化圖像總變差計(jì)算中用于各向異性加權(quán)的權(quán)重矩陣g，權(quán)重因子初始值均設(shè)為1。

（2）增加迭代次數(shù)k=k+ 1，利用代數(shù)迭代法重建得到更新的重建圖像u（k）。

（3）根據(jù)步驟（2）中得到的更新的重建圖像，求解滿足圖像總變差最小化的圖像。

（4）對(duì)選擇的圖像進(jìn)行各向異性的邊緣檢測，得到x和y 2個(gè)方向的邊緣圖像Ex以及Ey。

（5）根據(jù)邊緣圖像更新權(quán)重矩陣為gnew。

（6）使用更新后的權(quán)重矩陣計(jì)算圖像的圖像總變差。

重復(fù)步驟（2）～（6）直至滿足迭代截止條件，并輸出圖像。

圖1 Aniso-EGTV算法流程圖

圖2 Shepp-Logan理想圖像（true）

2 結(jié)果

對(duì)Shepp-Logan理想模型離散為256×256的矩陣（記為true，如圖2所示）后進(jìn)行等距扇形束投影，射線源到探測器的距離為800 mm，到旋轉(zhuǎn)中心的距離為400mm。探測器像素大小為0.806 6 mm，共有512個(gè)像素。360°掃描范圍內(nèi)均勻產(chǎn)生360幅投影，從投影中分別均勻抽取15、20、40個(gè)角度的投影數(shù)據(jù)，在無噪或添加高斯噪聲情況下進(jìn)行重建。

首先采用重建圖像做邊緣檢測，利用邊緣檢測結(jié)果做權(quán)重的更新，重建圖像如圖3的中間列所示，記作Aniso-EGTV-Mid。為了更好地得到邊緣信息，采用理想圖像做邊緣檢測，重建結(jié)果如圖3的最后一列所示，記作Aniso-EGTV-Prior。為了比較不同算法的重建效果，將采用TV-POCS方法重建的結(jié)果列于圖3第一列。同時(shí)，為了更好地觀察，放大虛線方框中感興趣區(qū)置于各圖的右下角。圖4畫出了各個(gè)稀疏角度下重建圖像第205行的剖線圖，可以看到Aniso-EGTV方法重建的結(jié)果與理想圖像更為接近，可以更好地恢復(fù)圖像的邊緣。

圖3 采用不同方法重建的圖像

3 討論

利用圖像梯度變換的稀疏化信息可以在較少的角度下重建CT圖像，且稀疏化模型直接影響著重建圖像的質(zhì)量。相比于傳統(tǒng)的TV-POCS方法，本文提出的Aniso-EGTV算法基于圖像各向異性的邊緣檢測信息更新TV計(jì)算中各分量的權(quán)重，用于后續(xù)的最優(yōu)化計(jì)算過程，從而實(shí)現(xiàn)各向異性的TV最小化重建方法。該方法可以更好地保持圖像邊緣的分辨率，重建結(jié)果與真實(shí)圖像更為接近。

圖4 采用不同方法重建的圖像第205行的剖線圖

Aniso-EGTV算法的重建效果依賴于邊緣檢測的精度，利用理想圖像邊緣檢測的結(jié)果重建時(shí)能夠減少由于數(shù)據(jù)不完全和噪聲導(dǎo)致的圖像變形。因此，好的邊緣檢測方法或者一些有關(guān)圖像邊緣的先驗(yàn)信息更有利于該方法的實(shí)施。但是好的邊緣檢測方法往往帶來參數(shù)選擇的困難，增加算法的復(fù)雜度。因此，實(shí)際中需要權(quán)衡邊緣檢測算法的實(shí)用性、復(fù)雜度、精度等之間的利害關(guān)系，選擇合適的邊緣檢測算法。

[1]Brenner D，Hall E.Computed tomography—an increasing source of radiation exposure[J].The New England Journal of Medicine，2007，

（????）（????）357（22）：2 277-2 284.

[2]Kalender W.Dose in x-ray computed tomography[J].Phys Med Biol，2014，59：R129-R150.

[3]Tuy H.An inverse formula for cone-beam reconstruction[J].SIAM J Appl Math，1983，43：546-552.

[4]Smith B.Image reconstruction from cone-beam projection：necessary and sufficient conditions and reconstruction methods[J].IEEE Trans Med Imaging，1985，4（1）：14-25.

[5]LIU Y，MA J，F(xiàn)AN Y，et al.Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction[J].Phys Med Biol，2012，57（23）：7 923-7 956.

[6]Sidky E Y，Kao C M，Pan X C.Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J].Journal of X-Ray Science and Technology，2006，14：119-139.

[7]Sidky E Y，Pan X.Image reconstruction in circular cone-beam computed tomography by constrained，total-variation minimization[J]. Phys Med Biol，2008，53：4 777-4 807.

[8]鄒晶，孫艷勤，張朋.由少量投影數(shù)據(jù)快速重建圖像的迭代算法[J].光學(xué)學(xué)報(bào)，2009，29（5）：1 198-2 004.

[9]Canny J.A computational approach to edge detection[J].IEEE Trans Pattern Anal Mach Intell，1986，8（6）：679-698.

（收稿：2014-03-26 修回：2014-06-17）

Anisotropic edge guided CT reconstruction from sparse views

RONG Jun-yan,LIU Wen-lei,GAO Peng,LIAO Qi-mei,LU Hong-bing
(School of Biomedical Engineering,the Fourth Military Medical University,Xi'an 710032,China)

ObjectiveTo develop a reconstruction algorithm for sparse view computed tomography (CT)to decrease the radiation dose.MethodsFirstly,the weights of the adding elements of the image total variation were updated according to the anisotropic edge detection information.Then the total variation minimization(TV minimization)algorithm was implemented to reconstruct the image.ResultsComputer simulation showed that the proposed algorithm could reconstruct the edge area much better than the conventional TV minimization.The precision of edge detection affected the reconstruction results.ConclusionThe proposed algorithm can enhance the image edge and improve the resolution and quality of reconstructed images.[Chinese Medical Equipment Journal，2014，35（10）：4-6，15]

anisotropic;edge detection;TV minimization;sparse view

R318；TH774

A

1003-8868（2014）10-0004-04

10.7687/J.ISSN1003-8868.2014.10.004

國家科技支撐計(jì)劃課題（2011BAI12B03）；中國博士后科學(xué)基金（2012M521863）；國家自然科學(xué)基金（81230035）

戎軍艷（1981—），女，博士，主要從事低劑量CT成像算法方面的研究工作，E-mail：junyanrong@126.com。

710032西安，第四軍醫(yī)大學(xué)生物醫(yī)學(xué)工程學(xué)院（戎軍艷，劉文磊，高 鵬，廖琪梅，盧虹冰）

盧虹冰，E-mail：luhb@fmmu.edu.cn