高大志，陳陽，楊剛，盧光明

1.南京軍區(qū)南京總醫(yī)院 醫(yī)學(xué)影像科，江蘇 南京 210002；2.東南大學(xué)生物醫(yī)學(xué)工程學(xué)院醫(yī)學(xué)信息研究所，江蘇 南京 210000

基于Bayesian非局部先驗(yàn)的低劑量X-射線CT成像

高大志1，陳陽2，楊剛1，盧光明1

1.南京軍區(qū)南京總醫(yī)院 醫(yī)學(xué)影像科，江蘇 南京 210002；2.東南大學(xué)生物醫(yī)學(xué)工程學(xué)院醫(yī)學(xué)信息研究所，江蘇 南京 210000

目的降低病人在CT掃描中所受的輻射劑量。方法在低劑量的掃描條件下，由于有效投影信息的缺少，CT圖像易受到量子噪聲的影響。本文提出了一種新的非局部先驗(yàn)Bayesian重建算法來提高CT圖像重建質(zhì)量。結(jié)果與傳統(tǒng)的圖像重建方法比較，本文提出的重建方法，能夠利用目標(biāo)圖像中更多形態(tài)結(jié)構(gòu)的全局信息來構(gòu)建先驗(yàn)項(xiàng)，從而克服了傳統(tǒng)圖像重建信息局限性的缺點(diǎn)。結(jié)論試驗(yàn)證明本文所提出的CT圖像重建方法在低劑量掃描條件下具有很好的表現(xiàn)。

螺旋CT；CT成像；Bayesian；低劑量掃描

本文導(dǎo)讀 >>

課題背景：目前低劑量X-射線CT掃描由于能夠減少病人的掃描輻射傷害而日益受到人們的關(guān)注。課題受國(guó)家自然科學(xué)基金（No.81000636）資助，主要針對(duì)低劑量CT掃描中由于有效光子數(shù)降低而導(dǎo)致成像質(zhì)量下降的問題所做的初步研究。

實(shí)驗(yàn)設(shè)計(jì)/論文構(gòu)思：基于Bayesian重建理論，通過引入圖像的全局信息來改善成像質(zhì)量，并與當(dāng)前國(guó)際上流行的其他圖像重建方法進(jìn)行比較和分析。

1 研究背景

X-射線CT掃描在臨床上的廣泛應(yīng)用，使得其中的輻射傷害日益受到人們的關(guān)注[1]。而目前螺旋CT掃描的普及也進(jìn)一步增加了病人在掃描中所受的輻射劑量[2]。因此，低劑量掃描條件(降低參數(shù)mAs、電壓kV)下的CT成像越來越越多地得到醫(yī)生和病人的重視[3]。然而，采用傳統(tǒng)的FBP重建算法，在低劑量掃描條件下, 投影數(shù)據(jù)g容易受到量子噪聲的影響，重建出的圖像f質(zhì)量會(huì)明顯下降，出現(xiàn)大量噪聲和星條狀偽影?；诮y(tǒng)計(jì)學(xué)的迭代重建能夠針對(duì)系統(tǒng)模型的物理效應(yīng)和探測(cè)數(shù)據(jù)、噪聲的統(tǒng)計(jì)泊松特性建立數(shù)學(xué)模型，同時(shí)可以通過Bayesian后驗(yàn)定理引入圖像空間的先驗(yàn)信息，解決圖像重建中由于噪聲而導(dǎo)致的病態(tài)問題，其重建的圖像質(zhì)量要優(yōu)于傳統(tǒng)的FBP方法[4-13]。Bayesian重建就是在迭代重建中引入圖像空間的先驗(yàn)信息項(xiàng)來抑制噪聲，基于此，我們可以構(gòu)建以下用于圖像重建的后驗(yàn)概率：

其中g(shù)和f分別為測(cè)量觀察數(shù)據(jù)和待重建的目標(biāo)衰減圖像。P(f)代表馬爾可夫隨機(jī)場(chǎng)(Markov Random Fields，MRF)先驗(yàn)分布。Z為正常數(shù)或者配分函數(shù)，為似然分布， 即為 MRF先驗(yàn)?zāi)芰炕蚰芰糠匠?，為像素點(diǎn)處的先驗(yàn)?zāi)芰糠匠蹋謪?shù)β控制MRF先驗(yàn)對(duì)重建圖像正則化作用的程度?？梢愿鶕?jù)(1)和(2)式建立相應(yīng)的后驗(yàn)?zāi)芰糠匠蹋?/p>

2005年，Buades等人設(shè)計(jì)了一種新的非局部去噪方法用于圖像的去噪處理[14]。在他們的非局部思想的啟發(fā)下，本文提出一個(gè)用于低劑量X-射線CT圖像Bayesian重建的非局部二次先驗(yàn)?zāi)Ｐ?，該模型不僅能使用目標(biāo)圖像中單個(gè)像素之間的灰度差信息，而且能夠有效的利用圖像中連接和連續(xù)的全局信息。相關(guān)的模擬試驗(yàn)表明：對(duì)于低劑量投影數(shù)據(jù)的X-射線CT圖像重建，本文所提出的非局部先驗(yàn)Bayesian重建法在降低噪聲效果和保持邊緣方面均具有非常好的表現(xiàn)。

1.1 傳統(tǒng)局部先驗(yàn)?zāi)Ｐ?/h3>

基于MRF理論，當(dāng)目標(biāo)圖像f滿足先驗(yàn)假設(shè)時(shí)，公式(2)中的能量方程U(f)具有最小值且相應(yīng)的先驗(yàn)分布(2)式達(dá)到最大值。而像素j點(diǎn)處的先驗(yàn)?zāi)芰糠匠蘒(f,j)通常等于以像素j與其鄰域Nj內(nèi)個(gè)點(diǎn)的差為自變量的勢(shì)能函數(shù)v( )的加權(quán)和：

可以通過選擇不同的勢(shì)能函數(shù)v(t)來設(shè)計(jì)用于Bayesian重建的不同先驗(yàn)?zāi)Ｐ?。如果?shì)能函數(shù)v(t)為二次的V(t)=t2形式，先驗(yàn)P（f）即為簡(jiǎn)單的二次QM平滑先驗(yàn)。我們也可通過使用選擇一些非二次勢(shì)能函數(shù)v(t)來使用非二次先驗(yàn)[3-10]，如Huber先驗(yàn)的勢(shì)能方程：

其中，δ為控制勢(shì)能函數(shù)的閾值參數(shù)。

公式(4)中的權(quán)值量wbj為表示圖像中像素b和j的相互關(guān)系的正常數(shù)。通常設(shè)定wbj的值同像素b和j之間距離成反比。下面即為廣泛應(yīng)用于圖像重建和恢復(fù)中的傳統(tǒng)局部先驗(yàn)的八鄰域歸一化權(quán)值圖和四鄰域歸一化權(quán)值圖:

以上的先驗(yàn)?zāi)Ｐ椭荒転閳D像重建提供固定的局部的先驗(yàn)信息。二次局部QM先驗(yàn)通過重建過程中的一個(gè)平均化的作用提供平滑的先驗(yàn)信息，從而易于導(dǎo)致圖像邊緣細(xì)節(jié)和噪聲同時(shí)平滑掉的過平滑效應(yīng)。具有邊緣保持作用的局部非二次先驗(yàn)則會(huì)因?yàn)闊o法有效的區(qū)分噪聲和一些較細(xì)微的細(xì)節(jié)而為圖像重建帶來負(fù)面的偽影效果。

1.2 非局部二次先驗(yàn)?zāi)Ｐ?/h3>

在建立非局部先驗(yàn)的過程中，一方面需要選擇一個(gè)較大的鄰域來包含圖像中的更多幾何形態(tài)信息，另一方面通過計(jì)算像素b的鄰域和j的鄰域的一個(gè)相似性測(cè)度，而不是用以上1.1節(jié)所述的像素間二維空間距離的簡(jiǎn)單反比例量來計(jì)算權(quán)值wbj。根據(jù)以上的非局部先驗(yàn)設(shè)計(jì)思路，我們可以根據(jù)(6)～(11)式為Bayesian重建設(shè)計(jì)新的非局部先驗(yàn)。

式中UNL(f)為非局部先驗(yàn)的能量函數(shù)， 這里Nj通常選擇一個(gè)較大尺寸的鄰域以引進(jìn)更多的結(jié)構(gòu)信息，wbj設(shè)為（1-0）二值函數(shù)，代表先驗(yàn)中像素b和j之間的權(quán)值，nb和nj為設(shè)定的以像素b和j為中心的方形鄰域， f(nb)和f(nj)為此兩個(gè)鄰域中所有像素灰度值數(shù)組，‖f（nb）-f（nj）‖代表此兩個(gè)像素灰度值數(shù)組的加權(quán)歐幾里得距離。

在選擇一個(gè)較大鄰域的基礎(chǔ)上，此非局部先驗(yàn)不僅比較圖像中兩個(gè)像素的灰度值，而且通過比較此兩個(gè)像素鄰域的相似性來計(jì)算勢(shì)能函數(shù)中的權(quán)值量。如圖1所示，在每組圖中，左邊的圖為原圖，右邊的圖描繪了左圖中心點(diǎn)鄰域中各點(diǎn)在非局部先驗(yàn)?zāi)Ｐ椭袡?quán)值的取值情況?？梢钥闯觯瑱?quán)值一般分布于比較相似的結(jié)構(gòu)處，對(duì)于兩個(gè)具有相似周圍結(jié)構(gòu)的像素點(diǎn)，此兩個(gè)像素點(diǎn)在先驗(yàn)中的權(quán)重較大，該非局部先驗(yàn)?zāi)Ｐ湍軌蚩紤]到圖像中的一些較大的幾何結(jié)構(gòu)形態(tài)的全局信息，能夠?qū)D像重建中的病態(tài)問題提供更為有效的正則化。

圖1 權(quán)值分布圖

在使用文中所述的非局部先驗(yàn)?zāi)Ｐ偷那闆r下，目標(biāo)圖(每組圖中左圖)中的中心點(diǎn)的鄰域權(quán)值分布如每組圖中的右圖所示，設(shè)定nb和nj的大小分別為31×31和7×7時(shí)，(a) 當(dāng)中心點(diǎn)位于一條豎直的邊緣上時(shí)，權(quán)值分布在此豎直邊緣線上；(b) 當(dāng)中心點(diǎn)位于傾斜的邊緣時(shí)，權(quán)值分布在傾斜的邊緣線上；(c) 當(dāng)中心點(diǎn)位于特定紋理結(jié)構(gòu)時(shí)，權(quán)值分布在特定的紋理結(jié)構(gòu)上。

2 實(shí)驗(yàn)與分析

2.1 CT迭代重建

在X-射線CT重建中，我們可以把測(cè)量投影數(shù)據(jù)理解為服從獨(dú)立泊松分布的隨機(jī)變量。第i對(duì)投影數(shù)值gi可以認(rèn)為是一個(gè)期望的泊松隨機(jī)數(shù)，這個(gè)是一個(gè)自變量為衰減圖像的方程[7-8]。似然方程PL(g/f)定義為在已知f的情況下，得到測(cè)量數(shù)據(jù)g的概率?；诖?，可以得到以下各式：

其中N是探測(cè)器對(duì)的數(shù)目，j是目標(biāo)圖像像素的數(shù)目，αi,j是在理想條件下圖像像素點(diǎn)j被探測(cè)器對(duì)i探測(cè)到的幾何概率。f2(g/f)對(duì)應(yīng)于(1)中P(g/f)，其對(duì)數(shù)似然能量方程為：

可以通過對(duì)后驗(yàn)?zāi)芰喀椎淖畲蠡\(yùn)算來重建出目標(biāo)圖像f：

然而，由(7)～(11)式可知，此非局部先驗(yàn)?zāi)芰恐械臋?quán)值項(xiàng)wbj的取值由目標(biāo)圖像f所決定，增加了直接求導(dǎo)計(jì)算的難度，實(shí)際試驗(yàn)我們采用如下步驟的重建算法：

（1）設(shè)定迭代初始項(xiàng)f°。

（2）權(quán)值項(xiàng)w的更新：在構(gòu)建此非局部先驗(yàn)?zāi)芰繒r(shí)，對(duì)于目標(biāo)圖像f中的每一個(gè)像素對(duì)(fb，fj)，用迭代過程中的當(dāng)前值的和公式(7)～(11)來計(jì)算wbj。

（3）目標(biāo)圖像f項(xiàng)的更新：

2.2 模擬試驗(yàn)

在實(shí)驗(yàn)中,我們使用單排探測(cè)器CT，X-射線源和探測(cè)器距離設(shè)為100cm，旋轉(zhuǎn)中心到探測(cè)器距離設(shè)為40cm，相鄰探測(cè)器距離為1mm，掃描系統(tǒng)為具有367個(gè)徑向取樣和360個(gè)角采樣的系統(tǒng)。圖2 (a)為試驗(yàn)中使用的256×256體模數(shù)據(jù)圖像，表示一個(gè)Shepp-Logan體模截面圖，圖像素點(diǎn)取值范圍為0～240；圖2(b)為對(duì)應(yīng)的正弦圖（sinogram）投影數(shù)據(jù)?；诖四M的兩個(gè)不同劑量(投影數(shù)據(jù)數(shù)值總和分別為1.2×108和1.5×108)的投影數(shù)據(jù)，投影數(shù)據(jù)中均加入了10%服從泊松分布的隨機(jī)噪聲。轉(zhuǎn)換概率矩陣A為一個(gè)平行帶狀積分幾何模型，通過(12)式中的模型進(jìn)行計(jì)算。

圖 2 體模圖像(a)與對(duì)應(yīng)的平行束投影正弦圖(b)

圖3和圖4分別顯示了使用以上的兩個(gè)不同劑量投影數(shù)據(jù)下的重建結(jié)果，使用Ramp濾波器的FBP重建(FBP1重建)，使用Hann濾波器的FBP重建(FBP2重建)和使用二次局部QM平滑先驗(yàn)，非二次局部Huber先驗(yàn)，MRP以及非局部先驗(yàn)的Bayesian重建圖像。對(duì)于后面4種Bayesian重建，選取第250次迭代的重建圖 (在實(shí)驗(yàn)中，當(dāng)?shù)降?50次時(shí)，重建圖像已趨于穩(wěn)定)。從圖3和圖4可以看到 探測(cè)數(shù)據(jù)在FBP重建圖像中易出現(xiàn)大量噪聲和星條偽影，相對(duì)于使用Ramp濾波器的FBP1重建(圖3 (a)和圖4 (a))；使用Hann濾波器的FBP2重建(圖3 (b) 和圖4(b))能夠更好的抑制噪聲；采用QM先驗(yàn)重建的圖像(圖3 (c) 和圖4(c))中出現(xiàn)了導(dǎo)致結(jié)構(gòu)信息模糊的過平滑效應(yīng)；對(duì)于使用非二次Huber先驗(yàn)的重建(圖3 (d) 、圖4 (d))和MRP(圖3 (e) 、圖4 (e))的重建來說，可以從結(jié)果圖像中觀察到塊狀的不規(guī)則偽影區(qū)域和階梯狀的偽影。另一方面，可以清楚的觀察到，使用本文提出的非局部先驗(yàn)重建的圖像(圖3 (f)和圖4 (f))具有更加均勻的背景區(qū)域和更加清晰的邊緣。采用該非局部先驗(yàn)的重建不僅能夠克服QM平滑先驗(yàn)的過平滑效應(yīng)，而且能夠在很大程度上解決非二次Huber先驗(yàn)和MRP所導(dǎo)致的階梯狀偽影的問題。

圖 3 使用總和為1.2×108探測(cè)數(shù)據(jù)sinogram的重建結(jié)果

圖 4 使用總和為1.5×108探測(cè)數(shù)據(jù)sinogram的重建結(jié)果

在重建試驗(yàn)中，使用不同的體模數(shù)據(jù)需要選擇不同的參數(shù)，試驗(yàn)中依據(jù)產(chǎn)生最高信噪比SNR（計(jì)算公式見(16)）的原則手工設(shè)定參數(shù)，試驗(yàn)中，信噪比SNR的計(jì)算公式為：

對(duì)于此兩種不同計(jì)數(shù)量的sinogram數(shù)據(jù)，在二次QM先驗(yàn)重建中，設(shè)定全局參數(shù)β的值分別為2×10-2和1.8×10-2；在采用非二次Huber先驗(yàn)的重建中，全局參數(shù)的值同QM一樣，設(shè)定全局參數(shù)δ的值分別為1.6和1.3；在采用MRP的重建中，設(shè)定全局參數(shù)β的值分別為40和20；在采用非局部先驗(yàn)的重建中，設(shè)定全局參數(shù)β的值分別為0.5和0.6，式(6)～(8)中參數(shù)n的取值分別為0.6和1.1，設(shè)定式(6)中nj的大小為11×11，設(shè)定(9)～(10)中的兩個(gè)比較鄰域nb和nj的大小為7×7。

表1給出了對(duì)于此兩種不同計(jì)數(shù)量的sinogram數(shù)據(jù)。由以上所有重建圖像對(duì)于真實(shí)體模圖像的信噪比SNR，可以看出使用本文提出的方法的重建圖像具有更高的信噪比。表2給出了不同重建方法在P4 1.6 GHz、512Mb RAM 的硬件環(huán)境下，matlab7.0的軟件環(huán)境下運(yùn)行得到重建結(jié)果所需的CPU時(shí)間?？梢钥闯霰疚乃岢龅姆蔷植肯闰?yàn)Bayesian重建需要更多的計(jì)算時(shí)間。

表1 兩種不同總探測(cè)量的sinogram數(shù)據(jù)，重建圖像相對(duì)于圖2中真實(shí)體模圖像數(shù)據(jù)的信噪比

表2 重建中使用不同方法所需CPU時(shí)間(s)

3 結(jié)論

相對(duì)于通常意義上的二次或者非二次的局部先驗(yàn)，本文提出的新的非局部先驗(yàn)?zāi)軌蚶媚繕?biāo)圖像中更多形態(tài)結(jié)構(gòu)的全局信息來構(gòu)建先驗(yàn)項(xiàng)，能夠引入更多更全面的先驗(yàn)信息，從而克服了局部先驗(yàn)信息局限性的缺點(diǎn)。基于模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)的試驗(yàn)表明：相對(duì)于使用普通的二次平滑先驗(yàn)及非二次邊緣保持先驗(yàn)的重建，使用該混合先驗(yàn)重建的CT圖像能夠更好保持邊緣信息，消除背景噪聲，且具有更高的信噪比。

進(jìn)一步的工作包括，對(duì)圖像重建中由非局部先驗(yàn)引進(jìn)的參數(shù)進(jìn)行有效的估計(jì)，減少使用非局部先驗(yàn)重建的計(jì)算量；應(yīng)使用更多的真實(shí)CT圖像掃描數(shù)據(jù)來測(cè)試該非局部先驗(yàn)重建方法在臨床試驗(yàn)中的效果。

[1]Linton OW, Mettler Jr FA. National conference on dose reduction in CT, with an emphasis on pediatric patients[J]. Am JRoentgenol,2003,181:321-329.

[2]Jung K, Lee K, Kim S, et al. Low-dose, volumetric helical CT:image quality, radiation dose, and usefulness for evaluation of bronchiectasis[J]. Invest Radiol,2000,35:557-563.

[3]Lu H, Hsiao I, Li X, et al. Noise properties of low-dose CT projections and noise treatment by scale transformations[C].2001 IEEE Nuclear Science Symp.Conf.,2001.p.1662-1666.

[4]Beekman FJ,Kamphuis C.Ordered subset reconstruction for X-ray CT[J]. PhysMed Biol,2001,46:1835-1855.

[5]Stan Z. Li. Markov Random Field Modeling in image Analysis[M]. Tokyo: Springer-Verlag, 2001:1-40.

[6]Black. M. J and Rangarajan A. Unification of line process,outlier rejection, and robust statistics with application in early vision [J]. Int. Journal of Computer Vision, 1996, 9:57-91.

[7]Nuyts J, De Man B, Dupont P, et al. Iterative reconstruction for helical CT:a simulation study[J].Phys Med Biol,1998,43:729-737.

[8]Sukovic P, Clinthorne NH. Penalized weighted leastsquares image reconstruction in single and dual energy X-ray computed tomography[J].IEEE Trans Med Imaging,2000,19(11):1075-1081.

[9]Lange K.Convergence of EM image reconstruction algorithms with Gibbs smoothness[J]. IEEE Trans Med Imaging,1990,9:439-446.

[10]Erdo-gan H, Fessler JA.Monotonic algorithms for transmission tomography[J]. IEEE Trans Med Imaging,1999,18(9):801-814.

[11]Alenius S, Ruotsalainen U, Astola J. Attenuation correction for PET using countlimited transmission images reconstructed with median root prior[J]. IEEE Trans Nucl Sci,1999,46:646-651.

[12]Yu DF, Fessler JA.Edge-preserving tomographic reconstruction with nonlocal regularization[J].IEEE Trans Med Imaging,2002,21(2):159-173.

[13]Li SZ. Markov random field modeling in image analysis[M].Tokyo: Springer-Verlag,2001:1-30.

[14]Buades A, Coll B, Morel JM. A nonlocal algorithm for image denoising[J]. Proc IEEE Int Conf Comput Vision Pattern Recognit,2005,2:60-65.

[15]Chen Y,Ma J,Feng Q,et al.Nonlocal prior Bayesian tomographic reconstruction[J]. J Math Imaging Vision,2008,30:133-146.

[16]周顯國(guó),等.貝葉斯決策分析在醫(yī)學(xué)步態(tài)分析中運(yùn)動(dòng)目標(biāo)檢測(cè)的應(yīng)用研究[J].中國(guó)醫(yī)療設(shè)備,2010,25(9):16-19.

Low-dose X-ray Computed Tomography Based on Bayesian Nonlocal Prior

GAO Da-zhi1, CHEN Yang2,YANG Gang1, LU Guang-ming1
1.Medical Imaging Department, Nanjing General Hospital of Nanjing Military Area Command, Nanjing Jiangsu 210002,China;2.Institute of Medical Information&Technology, School of Biomedical Engineering, Southeast University, Nanjing Jiangsu 210000,China

R319；TN919.81

A

10.3969/j.issn.1674-1633.2011.02.006

1674-1633(2011)02-0021-04

2010-10-20

2010-11-23

國(guó)家自然科學(xué)基金（81000636）資助。

本文作者：高大志，碩士研究生。

作者郵箱：njdazhi@hotmail.com

Abstract:ObjectiveTo reduce the radiation dose delivered to the patients during the application of computed tomography (CT).MethodsCT images can be easily degraded by the quantum noise under low X-ray dose scan protocols. Statistical reconstructions outperform the traditional filtered back-projection(FBP) reconstructions by accurately modeling the scan system and the measurement statistics. This paper aims to improve the CT reconstruction using a new nonlocal prior statistical reconstruction approach.ResultsCompared to traditional reconstruction approaches, the proposed nonlocal prior can impose an effective regularization for reconstructions by exploiting the image global information adaptively.ConclusionExperimentation validates that the proposed CT reconstructions have excellent performance with low-dose scan protocols.

Key words:spiral CT; CT imaging; Bayesian arithmetic; low-dose scanning