王中昊,夏 竟,李世杰,蔡志平

(國(guó)防科技大學(xué)計(jì)算機(jī)學(xué)院,湖南 長(zhǎng)沙 410073)

1 引言

隨著計(jì)算機(jī)斷層掃描CT(Computed Tomography)技術(shù)的不斷發(fā)展,錐形束斷層掃描CBCT (Cone Beam CT) 由于輻射劑量低、各項(xiàng)同性空間分辨率高等優(yōu)點(diǎn)[1],在臨床醫(yī)療、工業(yè)檢測(cè)等領(lǐng)域的應(yīng)用越來越廣泛。在不含金屬物體的情況下,CT裝置能夠拍攝出高質(zhì)量圖像,但如果有金屬假體等植入物出現(xiàn)在拍攝范圍FOV(Field Of View)內(nèi),將產(chǎn)生嚴(yán)重的金屬偽影,大幅降低圖像的質(zhì)量和臨床診斷價(jià)值。如何去除金屬偽影對(duì)圖像進(jìn)行恢復(fù),在CT研究中有著重要的現(xiàn)實(shí)意義。

目前,有許多針對(duì)扇形束CT的去金屬偽影算法MAR(Metal Artifact Reduction)。錐束CT的相關(guān)去金屬偽影算法大部分都在此基礎(chǔ)上提出的,原理上較為相似[2],即在原始圖像中將金屬區(qū)域分割后,采用一定方式對(duì)金屬區(qū)域進(jìn)行修復(fù),然后重建得到無金屬偽影重建圖像,再將原始圖像中的金屬部分置入,得到偽影校正圖像。這一過程中,最為關(guān)鍵的步驟是將金屬區(qū)域進(jìn)行準(zhǔn)確分割和對(duì)金屬區(qū)域進(jìn)行修復(fù)。應(yīng)用較為廣泛的算法是根據(jù)金屬與其他物質(zhì)的體素值的不同,設(shè)定閾值進(jìn)行分割,然后使用線性插值等插值算法對(duì)金屬區(qū)域進(jìn)行補(bǔ)全。雖然這些算法能夠在一定程度上抑制金屬偽影,但在實(shí)際應(yīng)用中效果仍不夠理想,甚至?xí)胄碌拇渭?jí)偽影,影響圖像質(zhì)量。

本文提出一種有效實(shí)用的算法,最大限度削減金屬偽影。該算法通過設(shè)計(jì)有效的金屬分割和補(bǔ)全方法,還原圖像真實(shí)結(jié)構(gòu),提升圖像質(zhì)量,同時(shí)減少二次偽影的引入。

本文的主要工作總結(jié)如下:

(1)提出了一種基于體素均值方差比VMR(Variance Mean Ratio)的金屬分割算法,通過體素點(diǎn)對(duì)應(yīng)投影值的變化特性對(duì)體素進(jìn)行區(qū)分,從而實(shí)現(xiàn)對(duì)金屬體素的精準(zhǔn)分割。

(2)在歸一化金屬偽影校正法NMAR(Normalized Metal Artifact Reduction)的基礎(chǔ)上,提出了一種基于雙調(diào)和方程插值修復(fù)的金屬偽影校正算法BIH-MAR(BIHarmonic Metal Artifact Reduction),能夠平滑地補(bǔ)全待修補(bǔ)區(qū)域,減少了次級(jí)偽影的引入,減小了對(duì)原始圖像細(xì)節(jié)的破壞,實(shí)現(xiàn)了對(duì)金屬偽影的消減和修復(fù)。

(3)在真實(shí)拍攝的CBCT圖像上進(jìn)行了實(shí)驗(yàn),并與常見的去金屬偽影算法進(jìn)行了對(duì)比,結(jié)果表明,本文提出的去金屬偽影算法效果優(yōu)于其他常用的算法,能夠較好地抑制金屬偽影,提升圖像質(zhì)量。

2 相關(guān)工作

在CBCT成像過程中,投影物體中如果包含金屬,則投影圖像在三維重建后,金屬區(qū)域周圍通常會(huì)產(chǎn)生嚴(yán)重的條紋狀痕跡,即金屬偽影。目前主流的去金屬偽影算法MAR主要分為4類[3]:投影插值法、迭代重建法、插值迭代混合法和深度學(xué)習(xí)校正法。

投影插值法算法因其直觀簡(jiǎn)單、計(jì)算量小而被廣泛應(yīng)用,其基本原理是對(duì)投影數(shù)據(jù)中的金屬部分,通過插值等方式進(jìn)行替換,重建后得到無金屬重建圖像,再與原始重建圖像中金屬部分相結(jié)合得到去偽影圖像。Lewitt等人[4]在1978年首先提出了投影插值法,采用多項(xiàng)式插值方法補(bǔ)全空心投影,同時(shí)對(duì)截?cái)嗤队斑M(jìn)行平滑連續(xù)性修復(fù)。Kalender等人[5]提出了針對(duì)金屬偽影的線性插值校正算法LIB-MAR(Linear Interpolation Based Metal Artifact Reduction)。該算法通過搜索算法在原始數(shù)據(jù)中分割出金屬跡線,再采用線性插值方法對(duì)金屬跡線區(qū)域進(jìn)行插值修補(bǔ),最后用反投影算法對(duì)修補(bǔ)后的正弦圖進(jìn)行重建。這種算法能夠較好地去除金屬偽影,是后續(xù)許多改進(jìn)算法的基礎(chǔ)算法,但該算法不可避免地會(huì)引入部分新的偽影,去偽影效果有限。Meyer等人[6]提出了一種歸一化去金屬偽影校正法NMAR,該算法基于先驗(yàn)圖像信息對(duì)原始數(shù)據(jù)進(jìn)行歸一化,從而使插值區(qū)域與周邊區(qū)域更加平滑,減少了次生偽影的出現(xiàn)。Liu等人[7]在NMAR的基礎(chǔ)上,針對(duì)CBCT設(shè)計(jì)出NMAR3,實(shí)現(xiàn)了較好的去偽影效果。Meilinger等人[8]在重建空間中對(duì)金屬及部分偽影進(jìn)行修正,然后對(duì)其進(jìn)行正向投影,再將正向投影與原始數(shù)據(jù)進(jìn)行融合得到去偽影圖像。這種算法更好地保留了原始投影中的細(xì)節(jié),但由于正向投影圖像與原始圖像存在一定截?cái)喱F(xiàn)象,融合過程中同樣會(huì)產(chǎn)生新的偽影。

迭代重建法首先對(duì)給定的一幅初始圖像進(jìn)行前向投影,然后將該投影與原始投影數(shù)據(jù)進(jìn)行比較,將其差值作為需要校正的誤差,不斷對(duì)圖像進(jìn)行校正。通過不斷迭代計(jì)算,對(duì)圖像信息不斷進(jìn)行檢驗(yàn)和修正直至誤差達(dá)到要求范圍。Gordon等人[9]首次將代數(shù)重建法ART(Algebraic Reconstruction Technique)引入圖像重建領(lǐng)域,其根據(jù)投影數(shù)據(jù)建立方程組,然后不斷迭代修正該方程組,從而得到重建區(qū)域的衰減系數(shù)。Levakhina等人[10]提出了一種基于加權(quán)迭代的校正算法,采用不同權(quán)重對(duì)投影數(shù)據(jù)進(jìn)行控制,并根據(jù)投影之間的差異度設(shè)置加權(quán)因子,使用加權(quán)因子控制反投影數(shù)值。當(dāng)射線穿過金屬時(shí),其本身差異度較高,加權(quán)因子設(shè)置較低,在迭代中貢獻(xiàn)也較低,從而抑制了金屬偽影的產(chǎn)生。這些迭代重建算法能夠較好地完成對(duì)金屬偽影的校正,然而由于其計(jì)算步驟繁瑣,所需時(shí)間較長(zhǎng)等問題,在實(shí)際應(yīng)用中的執(zhí)行難度較大。

Figure 1 Flow chart of metal artifact correction in cone beam CT圖1 錐束CT金屬偽影校正流程圖

隨著計(jì)算能力的提升和深度學(xué)習(xí)技術(shù)的發(fā)展,金屬偽影校正領(lǐng)域近年來也出現(xiàn)了許多基于深度學(xué)習(xí)的去偽影算法。Gjesteby等人[11]首次提出將卷積神經(jīng)網(wǎng)絡(luò)CNN(Convolutional Neural Network)和NMAR算法進(jìn)行融合,從而進(jìn)一步校正NMAR處理后的金屬偽影部分,進(jìn)一步改善了NMAR的校正效果。Ghani等人[12]提出了一種基于生成對(duì)抗網(wǎng)絡(luò)GAN(Generative Adversarial Networks)的金屬偽影校正算法,通過訓(xùn)練網(wǎng)絡(luò)對(duì)金屬部分的投影數(shù)據(jù)進(jìn)行補(bǔ)全,再通過濾波反投影[13]算法進(jìn)行重建?？梢钥吹?基于深度學(xué)習(xí)的算法在特定場(chǎng)景下能夠取得較好的效果,但由于偽影的形成根據(jù)不同的拍攝條件和金屬物體的性質(zhì)具有較大不確定性,導(dǎo)致其泛化能力不佳且訓(xùn)練樣本不易獲得,同樣存在一定局限性。

本文在歸一化去金屬偽影校正法的基礎(chǔ)上,結(jié)合基于方差均值比的金屬閾值分割法和基于雙調(diào)與方程插值的圖像修復(fù)算法,實(shí)現(xiàn)錐形束CT的金屬偽影校正。相比現(xiàn)在已有的算法,本文所提算法能夠在滿足實(shí)用性的基礎(chǔ)上,提升金屬偽影校正效果。

3 本文算法

本文提出的錐束CT金屬偽影校正算法流程如圖1所示。

首先,對(duì)含金屬偽影的重建圖像進(jìn)行雙邊濾波,并對(duì)金屬區(qū)域進(jìn)行分割,獲得金屬圖像和去金屬圖像。然后,對(duì)金屬圖像和去金屬圖像進(jìn)行正向投影和歸一化操作,獲得金屬投影區(qū)域和歸一化投影。接下來,對(duì)歸一化投影中的金屬投影區(qū)域進(jìn)行雙調(diào)和插值修復(fù),得到修復(fù)的無金屬圖像。最后,對(duì)修復(fù)的金屬圖像進(jìn)行去歸一化和FDK(Feldkamp-Davis-Kress)[14,15]算法重建,得到無金屬重建圖像,并與原始金屬圖像進(jìn)行融合,獲得最終的校正圖像。

3.1 圖像雙邊濾波

原始圖像經(jīng)過FDK重建后,仍存在一定噪聲。為了更好地去除噪聲同時(shí)保留圖像邊緣結(jié)構(gòu),采用雙邊濾波器BF(Bilateral Filter)[16]對(duì)圖像進(jìn)行濾波,以有效地去除噪聲和平滑圖像,并保留圖像的邊緣和細(xì)節(jié)特征。BF是一種基于空間域和灰度值域的濾波器。對(duì)于一個(gè)中心像素,BF考慮其周圍像素的空間距離和像素值之間的相似性,用高斯權(quán)值函數(shù)進(jìn)行加權(quán)平均。其中,空間高斯權(quán)值函數(shù)根據(jù)像素之間的空間距離進(jìn)行計(jì)算,像素高斯權(quán)值函數(shù)則根據(jù)像素之間的差異進(jìn)行計(jì)算。通過這種方式,雙邊濾波器可以在保留圖像細(xì)節(jié)的同時(shí)去除噪聲。設(shè)V′(x,y,z)為BF濾波后的體素值,其計(jì)算如式(1)所示:

(1)

其中,V(x,y,z)為原始重建圖像中位置(i,j,k)處的像素值;Np是以像素點(diǎn)(x,y,z)為中心的鄰域像素集合;w(i,j,k,x,y,z)是雙邊濾波器中的權(quán)重函數(shù),用來衡量像素(i,j,k)與像素(x,y,z)之間的相似度,其通常包括一個(gè)空間權(quán)重和一個(gè)像素權(quán)重,如式(2)所示:

w(i,j,k,x,y,z)=

wspace(i,j,k,x,y,z)×wpixel(i,j,k,x,y,z)

(2)

其中,空間權(quán)重考慮了像素之間的空間距離,其定義如式(3)所示:

wspace(i,j,k,x,y,z)=

(3)

其中,σd是控制空間權(quán)重衰減速度的參數(shù),本文根據(jù)實(shí)踐結(jié)果設(shè)置為濾波器卷積核半徑。像素權(quán)重考慮了像素之間的相似度,采用高斯函數(shù)來定義,如式(4)所示:

wpixel(i,j,k,x,y,z)=

(4)

其中,σr是控制像素權(quán)重衰減速度的參數(shù),本文設(shè)置為一個(gè)較小的值。通過控制像素間的空間距離和灰度變化范圍調(diào)節(jié)像素的加權(quán)值,實(shí)現(xiàn)對(duì)圖像的濾波。

3.2 基于體素方差均值比的金屬分割

在CT影像中,使用CT值衡量組織密度,使用亨氏單位HU (Hounsfield Unit)作為CT值的計(jì)量單位。在重建圖像中,準(zhǔn)確對(duì)金屬區(qū)域分割是金屬偽影校正算法的基礎(chǔ)。基于閾值的金屬分割算法具有計(jì)算簡(jiǎn)單、性能穩(wěn)定等優(yōu)點(diǎn),應(yīng)用較為廣泛。然而,金屬與其周圍的非金屬物體的體素變化在重建體積中是連續(xù)的,位于金屬物體邊緣的非金屬物體受偽影影響,其CT值接近金屬的CT值,難以設(shè)定一個(gè)特定的CT值將金屬進(jìn)行準(zhǔn)確地分割。例如,一般的金屬CT值遠(yuǎn)大于2 000 HU,當(dāng)閾值設(shè)定為2 000 HU時(shí),金屬周圍部分非金屬體素在偽影的影響下也會(huì)超過閾值從而被分割為金屬;當(dāng)閾值設(shè)定為4 000 HU時(shí),部分金屬邊緣又會(huì)因其體素小于閾值而未被正確分割。因此,為了更加精準(zhǔn)地對(duì)金屬進(jìn)行分割,本文引入體素對(duì)應(yīng)投影點(diǎn)的方差均值比VMR(Variance-to-Mean Ratio)作為判斷依據(jù),輔助分割。

方差均值比即變異系數(shù),是衡量觀測(cè)值變異程度的一個(gè)統(tǒng)計(jì)量。本文將該比率用于CBCT 重建過程中,作為某一體素受到CT偽影影響概率的度量。方差均值比定義為某一體素在所有角度下原始投影圖像對(duì)應(yīng)的投影點(diǎn)值的方差除以該體素的CT值。設(shè)給定重建圖像位置i處的體素CT值為V(i),該體素的中心點(diǎn)被射線源投影到投影圖像pθ的位置qi,θ處,其中θ∈(0,π)為旋轉(zhuǎn)角度。則體素點(diǎn)i的方差均值比VMR(i)定義如式(5)所示:

(5)

(6)

其中,Tmetal為金屬閾值,Tvmr為VMR閾值,可通過金屬物占比結(jié)合直方圖法確定其取值。當(dāng)某一體素點(diǎn)的CT值大于或等于閾值Tmetal,且VMR值小于閾值Tvmr時(shí),該體素點(diǎn)判定為金屬并將其在二值圖像中的像素值設(shè)為1,其他區(qū)域的像素值設(shè)為0,然后對(duì)金屬區(qū)域進(jìn)行前向投影,得到金屬物體在投影平面中的區(qū)域。

3.3 生成歸一化圖像

直接對(duì)金屬區(qū)域進(jìn)行插值修補(bǔ)會(huì)因?yàn)榇扪a(bǔ)區(qū)域周圍不平滑而引入次生偽影[17],因此本文采用歸一化金屬偽影校正,需要對(duì)原始重建圖像進(jìn)行前向投影生成先驗(yàn)圖像。為了減少金屬偽影對(duì)非金屬物的影響,首先要將原始重建圖像進(jìn)行量化,將金屬區(qū)域填充為空氣的CT值,其他區(qū)域分別填充為對(duì)應(yīng)物體的平均CT值。本文實(shí)驗(yàn)中使用陶瓷杯和金屬螺釘進(jìn)行實(shí)驗(yàn),可將原始重建圖像量化為金屬、陶瓷和空氣區(qū)域,并對(duì)不同物質(zhì)進(jìn)行填充賦值,然后將其映射為對(duì)應(yīng)的二元CT值,如式(7)所示:

(7)

(8)

3.4 雙調(diào)和方程插值補(bǔ)全

Δ2P(x,y)=u(x,y)

(9)

其中,Δ2表示拉普拉斯算子的平方,對(duì)歸一化圖像中待修補(bǔ)的金屬區(qū)域使用雙調(diào)和插值;P(x,y)表示待修補(bǔ)的圖像區(qū)域;u(x,y)表示已知的金屬區(qū)域邊緣。該方程是要尋找一個(gè)能夠充分平滑并且與已知邊緣相符合的函數(shù)P(x,y)。

采用Canny算子從原始圖像中提取出邊緣信息后,將雙調(diào)和方程代入求解器中,迭代求解得到待修復(fù)區(qū)域P(x,y)。雙調(diào)和方程的求解公式如式(10)所示:

(10)

其中,G(x,y,s,t)是點(diǎn)源函數(shù),表示在點(diǎn)(s,t)處放置一個(gè)單位源,產(chǎn)生的響應(yīng)在點(diǎn)(x,y)處的值是該偏微分方程在特定邊界條件下的響應(yīng);u(s,t)是已知的邊緣信息;log(1/r)是平滑參數(shù),r取值較小的圖像就會(huì)更接近原圖像,但可能會(huì)存在較多的噪聲和細(xì)節(jié),反之圖像會(huì)更加平滑,但可能導(dǎo)致邊緣信息丟失。實(shí)踐證明,可以適當(dāng)增大r的取值,以保證平滑度,減少次級(jí)偽影。

Figure 2 Slices of original reconstructed image圖2 原始重建圖像切片

(11)

最后,使用FDK算法對(duì)P′θ進(jìn)行重建,得到不含金屬的重建圖像,將原始重建圖像中的金屬部分置入,得到最終的偽影消減圖像,其融合過程如式(12)所示:

(12)

4 實(shí)驗(yàn)與結(jié)果分析

實(shí)驗(yàn)數(shù)據(jù)通過實(shí)驗(yàn)機(jī)器實(shí)拍獲得,采用CBCT掃描方式,通過機(jī)架繞其水平軸旋轉(zhuǎn)對(duì)實(shí)驗(yàn)物體進(jìn)行正向投影。實(shí)驗(yàn)使用的射線源到平板探測(cè)器中心距離為600 mm,到旋轉(zhuǎn)中心的距離為400 mm;射線源管電壓為110 kV,電流強(qiáng)度為10.9 mA;平板探測(cè)器的像素矩陣尺寸為1274×1024,像素大小為0.125 mm×0.125 mm,掃描范圍為0～360°,步長(zhǎng)為0.5°,機(jī)架旋轉(zhuǎn)掃描一周得到原始投影圖像720幅。實(shí)驗(yàn)選取陶瓷杯和金屬螺釘進(jìn)行掃描,首先對(duì)陶瓷杯進(jìn)行單獨(dú)掃描作為對(duì)照數(shù)據(jù),然后將金屬螺釘附著在陶瓷杯上進(jìn)行掃描得到實(shí)驗(yàn)數(shù)據(jù)。

采用FDK算法對(duì)采集的投影數(shù)據(jù)(圖2a)進(jìn)行三維圖像重建(圖2b和圖2c)。從圖2可以看出,重建圖像中可見明顯的呈明亮條形和暗帶的金屬偽影,圖像質(zhì)量較差。

Figure 3 Images at each stage of the proposed correction algorithm圖3 本文校正算法各階段圖像

采用本文提出的BIH-MAR算法進(jìn)行校正。圖3給出了算法的各階段圖像:(1)對(duì)重建圖像進(jìn)行雙邊濾波去除部分圖像噪聲,圖3a顯示了濾波處理后的圖像;(2)結(jié)合閾值和VMR對(duì)圖像金屬部分進(jìn)行分割,圖3b為進(jìn)行金屬分割后得到的二值化圖像對(duì)比圖,上半部分為原始切片,下半部分為金屬部分的二值分割,可以看到金屬分割準(zhǔn)確清晰;(3)將原始重建圖像化為二值CT后進(jìn)行前向投影,得到先驗(yàn)圖像如圖3c所示;(4)使用先驗(yàn)圖像對(duì)原始投影圖像進(jìn)行歸一化處理,對(duì)該歸一化圖像的金屬部分采用雙調(diào)和方程進(jìn)行插值補(bǔ)全,得到修復(fù)后的投影圖像如圖3d所示;(5)使用修復(fù)的投影圖像進(jìn)行FDK重建,得到無金屬偽影校正圖像如圖3e所示;(6)將原始重建圖像的金屬部分與無金屬偽影校正圖像進(jìn)行融合,得到最終的校正圖像如圖3f所示。分別采用基于二次線性插值的去金屬偽影算法LIB-MAR[5]和歸一化去金屬偽影算法NMAR[6]對(duì)原始重建圖像進(jìn)行校正,與本文BIH-MAR算法的結(jié)果進(jìn)行對(duì)比,如圖4所示。

Figure 4 Effect comparison of three metal artifact correction algorithms圖4 3種金屬偽影校正算法效果對(duì)比

從圖4可以看出,與原始重建圖像(圖2c)相比,LIB-MAR算法在消減了一定金屬偽影的同時(shí),也引入了大量的次級(jí)偽影,導(dǎo)致部分結(jié)構(gòu)被破壞,金屬偽影校正效果較差;NMAR算法效果有明顯提升,使金屬偽影得到了較好的抑制,原始結(jié)構(gòu)得到保護(hù),圖像質(zhì)量有一定提升,但仍存在一些偽影未得到抑制;與這2種算法相比,BIH-MAR算法能夠更大程度地對(duì)金屬偽影進(jìn)行消減,同時(shí)組織結(jié)構(gòu)保留完整,整體圖像質(zhì)量提升效果明顯,校正效果最優(yōu)。

為了對(duì)這3種MAR算法進(jìn)行定量比較[20],在偽影較為嚴(yán)重的Z方向第200層切片上選取3個(gè)尺寸為50×50的感興趣區(qū)域ROI(Region of Interest)(如圖5所示),以無金屬投影生成的重建圖像作為參考,計(jì)算ROI區(qū)域的均方根誤差RMSE(Root Mean Squared Error)進(jìn)行評(píng)估。RMSE的計(jì)算如式(13)所示:

(13)

Figure 5 Three ROI regions of slice of layer 200 in the Z direction圖5 Z向第200層切片的3個(gè)ROI區(qū)域

表1顯示了RMSE結(jié)果,在同一ROI上和所有ROI上,BIH-MAR算法對(duì)應(yīng)的RMSE值相對(duì)于其他2種算法的均為最小;在所有ROI上,BIH-MAR算法對(duì)應(yīng)的RMSE為0.028,比LIB-MAR和NMAR 算法的分別減少了8%和22%,表明BIH-MAR 算法對(duì)應(yīng)的校正圖像與原始圖像的偏差最小,對(duì)金屬偽影的校正效果最好,同時(shí)更好地保留了圖像原始結(jié)構(gòu),與其他2種算法相比更加有效。

Table1 RMSE values of ROI regions of different MAR algorithms表1 不同MAR算法ROI區(qū)域的RMSE值

5 結(jié)束語(yǔ)

本文提出了一種歸一化校正的自適應(yīng)錐束CT金屬偽影消減算法。該算法首先對(duì)原始重建圖像進(jìn)行雙邊濾波,去除圖像底噪,然后結(jié)合VMR對(duì)金屬進(jìn)行閾值分割,將原始圖像映射為二元CT后進(jìn)行前向投影,得到不含金屬的先驗(yàn)投影圖像;然后使用先驗(yàn)投影圖像對(duì)原始投影進(jìn)行歸一化,并對(duì)金屬區(qū)域采用雙調(diào)和方程插值補(bǔ)全,得到修復(fù)后的無金屬投影圖像;最后使用FDK算法對(duì)修復(fù)后的無金屬投影圖像進(jìn)行三維重建,并與金屬圖像融合,獲得最終的金屬偽影校正圖像。利用實(shí)拍數(shù)據(jù)對(duì)金屬偽影校正算法進(jìn)行有效性驗(yàn)證。實(shí)驗(yàn)結(jié)果表明,本文算法的金屬偽影校正效果優(yōu)于常用的LIB-MAR和NMAR算法的,從定性和定量的角度驗(yàn)證了本文算法的有效性。