基于雙調和插值的錐束CT金屬偽影校正算法*

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

(國防科技大學計算機學院,湖南 長沙 410073)

1 引言

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

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

本文提出一種有效實用的算法,最大限度削減金屬偽影。該算法通過設計有效的金屬分割和補全方法,還原圖像真實結構,提升圖像質量,同時減少二次偽影的引入。

本文的主要工作總結如下:

(1)提出了一種基于體素均值方差比VMR(Variance Mean Ratio)的金屬分割算法,通過體素點對應投影值的變化特性對體素進行區分,從而實現對金屬體素的精準分割。

(2)在歸一化金屬偽影校正法NMAR(Normalized Metal Artifact Reduction)的基礎上,提出了一種基于雙調和方程插值修復的金屬偽影校正算法BIH-MAR(BIHarmonic Metal Artifact Reduction),能夠平滑地補全待修補區域,減少了次級偽影的引入,減小了對原始圖像細節的破壞,實現了對金屬偽影的消減和修復。

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

2 相關工作

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

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

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

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

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

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

3 本文算法

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

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

3.1 圖像雙邊濾波

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

(1)

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

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

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

(2)

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

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

(3)

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

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

(4)

其中,σr是控制像素權重衰減速度的參數,本文設置為一個較小的值。通過控制像素間的空間距離和灰度變化范圍調節像素的加權值,實現對圖像的濾波。

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

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

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

(5)

(6)

其中,Tmetal為金屬閾值,Tvmr為VMR閾值,可通過金屬物占比結合直方圖法確定其取值。當某一體素點的CT值大于或等于閾值Tmetal,且VMR值小于閾值Tvmr時,該體素點判定為金屬并將其在二值圖像中的像素值設為1,其他區域的像素值設為0,然后對金屬區域進行前向投影,得到金屬物體在投影平面中的區域。

3.3 生成歸一化圖像

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

(7)

(8)

3.4 雙調和方程插值補全

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

(9)

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

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

(10)

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

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

(11)

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

(12)

4 實驗與結果分析

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

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

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

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

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

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

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

(13)

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

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

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

5 結束語

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