劉仰川，高 鵬,2，朱葉晨，高 欣

一種基于先驗(yàn)圖像的錐束CT金屬偽影校正算法

劉仰川1，高 鵬1,2，朱葉晨1，高 欣1

(1. 中國(guó)科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所醫(yī)學(xué)影像技術(shù)研究室，江蘇 蘇州 215163；2. 南京理工大學(xué)電子工程與光電技術(shù)學(xué)院，江蘇 南京 210094)

為了有效抑制錐束CT(CBCT)重建中金屬植入物引入的偽影，提出一種基于先驗(yàn)圖像的金屬偽影校正算法。首先對(duì)含金屬偽影的重建圖像進(jìn)行雙邊濾波、金屬閾值分割、組織聚類等預(yù)處理，獲得金屬圖像和不含金屬信息的先驗(yàn)圖像；再對(duì)二者正向投影，獲得金屬投影區(qū)域和先驗(yàn)投影數(shù)據(jù)；而后利用先驗(yàn)投影數(shù)據(jù)及金屬邊界鄰域的投影數(shù)據(jù)對(duì)金屬投影區(qū)域插值，獲得修復(fù)的投影數(shù)據(jù)；最后利用FDK算法對(duì)修復(fù)的投影數(shù)據(jù)重建，并將其與金屬圖像融合，獲得最終的校正圖像。為了驗(yàn)證該算法的性能，利用三維Shepp-Logan頭部模型數(shù)據(jù)和臨床頭部CT數(shù)據(jù)開展金屬偽影校正實(shí)驗(yàn)，結(jié)果表明：與常用的線性插值算法和圖像修補(bǔ)算法相比，該算法的校正圖像均方根誤差最小、峰值信噪比最大。這說明該算法在有效保留圖像邊緣信息的同時(shí)，可有效地抑制金屬偽影。

錐束CT；金屬偽影校正；雙邊濾波；先驗(yàn)圖像；插值

近年來，錐束CT(cone beam computed tomography，CBCT)因體積小、輻射劑量低、各向同性空間分辨率高等優(yōu)點(diǎn)[1]，在口腔科、骨科的應(yīng)用越來越廣泛，是目前最有前途和實(shí)用性的成像方式之一[2-3]。然而，病人體內(nèi)攜帶的金屬植入物[3](如金屬義齒、骨釘?shù)?會(huì)給重建圖像帶來嚴(yán)重的條狀偽影(金屬偽影)，極大影響診斷準(zhǔn)確性。金屬偽影校正[4](metal artifact reduction，MAR)算法可以在圖像重建過程中，減少甚至消除條狀偽影，改善重建圖像質(zhì)量。

目前，錐束CT的MAR算法主要分為3類[5]：①插值法：從原重建圖像或投影數(shù)據(jù)中分割出金屬區(qū)域，并對(duì)金屬區(qū)域插值，再將重建圖像與金屬圖像融合，實(shí)現(xiàn)金屬偽影校正；②迭代法：借助迭代重建算法的平滑和去偽影特性，采用不同參數(shù)對(duì)金屬投影區(qū)域和非金屬投影區(qū)域分別重建，再進(jìn)行圖像融合，實(shí)現(xiàn)金屬偽影校正；③混合法：將2種或多種MAR方法結(jié)合(通常是將插值法、迭代法結(jié)合)，以提升校正性能。迭代法和混合法復(fù)雜度高、耗時(shí)嚴(yán)重[6-7]，臨床應(yīng)用局限性大，而插值法復(fù)雜度低、計(jì)算效率高，是當(dāng)前MAR主流算法。因此，本文圍繞插值法開展研究。MEILINGER等[8]提出了一種基于互信息和邊緣濾波的MAR算法，利用閾值法從原重建圖像中分割出金屬圖像，并將金屬鄰域像素點(diǎn)插值到金屬區(qū)域，再通過正向投影及邊緣濾波，獲得修復(fù)的投影數(shù)據(jù)。ZHANG等[9]利用錐束CT的成像幾何關(guān)系定位金屬投影區(qū)域，再利用拉普拉斯擴(kuò)散法將金屬邊緣投影數(shù)據(jù)插值到金屬區(qū)域，獲得修復(fù)的投影數(shù)據(jù)。BEISTER等[10]提出了一種基于前投影數(shù)據(jù)的簡(jiǎn)約MAR算法，利用閾值法從原重建圖像中分割出金屬圖像，并對(duì)金屬圖像正向投影確定其投影位置，再利用三維線性插值法將金屬邊緣投影數(shù)據(jù)插值到金屬投影區(qū)域，獲得修復(fù)的投影數(shù)據(jù)。以上算法簡(jiǎn)單、快速，但校正圖像中容易出現(xiàn)嚴(yán)重的次級(jí)偽影(因校正不當(dāng)而產(chǎn)生的新偽影)。為此，PRELL等[11]提出了一種可降低次級(jí)偽影的MAR算法，先對(duì)原重建圖像高斯濾波和金屬閾值分割，再對(duì)分割出的金屬圖像正向投影，后使用三維線性插值算法對(duì)金屬投影區(qū)域插值，并進(jìn)行歸一化處理，獲得修復(fù)的投影數(shù)據(jù)。該算法能較好地消除金屬偽影并產(chǎn)生較少的次級(jí)偽影，但校正圖像中的組織邊緣容易被模糊，導(dǎo)致邊緣信息缺失。

針對(duì)現(xiàn)有插值類MAR算法校正后的圖像存在次級(jí)偽影及邊緣模糊的問題，本文提出一種基于先驗(yàn)圖像的CBCT金屬偽影校正算法(prior-image- based MAR，PIB-MAR)。該算法在2個(gè)方面做出改進(jìn)：①采用雙邊濾波對(duì)原重建圖像進(jìn)行預(yù)處理，去除噪聲并保留邊緣信息；②采用先驗(yàn)投影數(shù)據(jù)及金屬鄰域投影數(shù)據(jù)對(duì)金屬投影區(qū)域插值，抑制次級(jí)偽影產(chǎn)生。為了驗(yàn)證算法的性能，利用三維Shepp-Logan頭部模型數(shù)據(jù)和臨床頭部CT數(shù)據(jù)進(jìn)行金屬偽影校正實(shí)驗(yàn)，并與常規(guī)的線性插值算法(linear-interpolation-based MAR，LIB-MAR)[11]、圖像修復(fù)算法(image-inpainting-based MAR，IIB- MAR)[12-13]進(jìn)行比較。

1 錐束CT金屬偽影校正算法

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

第1步：預(yù)處理。對(duì)含金屬偽影的重建圖像進(jìn)行雙邊濾波、金屬閾值分割、組織聚類等預(yù)處理，獲得金屬圖像和先驗(yàn)圖像。

第2歩：金屬偽影校正。對(duì)金屬圖像和先驗(yàn)圖像進(jìn)行正向投影，獲得金屬投影區(qū)域和先驗(yàn)投影數(shù)據(jù)，再對(duì)金屬投影區(qū)域進(jìn)行修復(fù)，獲得修復(fù)的投影數(shù)據(jù)。

圖1 錐束CT金屬偽影校正流程圖

第3步：CT圖像重建。利用FDK算法對(duì)修復(fù)的投影數(shù)據(jù)進(jìn)行重建，得到中間重建圖像，與第1步分割出的金屬圖像進(jìn)行融合，獲得最終的校正圖像。

1.1 雙邊濾波

針對(duì)含金屬偽影的CT圖像，采用雙邊濾波器(bilateral filter，BF)[14-15]進(jìn)行濾波，去除部分偽影和噪聲。BF是一種非線性濾波器，基于像素間的歐式距離及灰度差異，在有效去除圖像噪聲的同時(shí)，保留圖像邊緣信息。針對(duì)CT圖像，BF采用高斯分布的加權(quán)平均函數(shù)對(duì)圖像進(jìn)行平滑，權(quán)重系數(shù)由2部分組成：①像素間的灰度差值范圍，即像素范圍域?yàn)V波核函數(shù)；②像素間的歐式距離，即空間域?yàn)V波核函數(shù)。

設(shè)()為含金屬偽影的重建圖像，將灰度值進(jìn)行歸一化，經(jīng)雙邊濾波，得到輸出圖像，即

1.2 金屬閾值分割

針對(duì)雙邊濾波后的CT圖像，采用閾值法進(jìn)行金屬分割，獲得金屬圖像。閾值分割是一種傳統(tǒng)且最常用的圖像分割算法[17-18]，具有簡(jiǎn)單、計(jì)算量小、性能穩(wěn)定等優(yōu)點(diǎn)。在MAR算法中，金屬分割的準(zhǔn)確性至關(guān)重要，可直接影響其效果。

CT圖像中，不同組織的CT值差異較大，如空氣的CT值為–1 000 HU、脂肪為–120~–90 HU、骨為300~2 000 HU，而各類金屬的CT值遠(yuǎn)大于2 000 HU，甚至上萬。因此，可采用閾值法將金屬從CT圖像中分割出來，即

其中，為閾值，可通過直方圖法確定取值，一般設(shè)為最大像素值的30%。從上式可以看出，像素值大于和等于的區(qū)域?yàn)榻饘賵D像，在二值圖像中設(shè)為1，其他區(qū)域的像素值設(shè)為0。

1.3 先驗(yàn)圖像生成

利用軟組織CT值對(duì)CT圖像中分割出的金屬區(qū)域進(jìn)行填充；借助“類組織模型”的思想[19]，采用三維K-means算法[20]將填充后的CT圖像聚類為空氣、脂肪、軟組織和骨；給不同聚類的組織賦值，獲得不含金屬信息的先驗(yàn)圖像prior，即

1.4 投影修復(fù)

根據(jù)錐束CT成像幾何關(guān)系，分別對(duì)分割出的金屬圖像、生成的先驗(yàn)圖像進(jìn)行正向投影，獲得金屬投影區(qū)域和先驗(yàn)投影數(shù)據(jù)。如果直接利用先驗(yàn)投影數(shù)據(jù)替換原投影數(shù)據(jù)中的金屬區(qū)域，修復(fù)后的邊緣連續(xù)性差，容易導(dǎo)致重建圖像中出現(xiàn)次級(jí)偽影。為此，采用改進(jìn)的修復(fù)方法，即利用先驗(yàn)投影數(shù)據(jù)、金屬邊界鄰域的投影數(shù)據(jù)對(duì)金屬投影區(qū)域進(jìn)行插值。

1.5 FDK重建與融合

采用FDK算法[21]對(duì)修復(fù)后的投影數(shù)據(jù)重建，獲得中間重建圖像，并與分割出的金屬圖像融合，獲得最終的校正圖像。FDK是在二維濾波反投影(filtered back-projection，F(xiàn)BP)算法基礎(chǔ)上發(fā)展出的一種三維重建算法。FDK屬于近似重建算法，其原理是將錐形束射線視為具有不同傾斜角度的扇形束射線面的集合，中心扇形束射線面(垂直于平板探測(cè)器)對(duì)應(yīng)的投影數(shù)據(jù)采用FBP算法重建，而非中心扇形束射線面對(duì)應(yīng)的投影數(shù)據(jù)采用修正的FBP算法重建，從而將錐形束重建問題轉(zhuǎn)化為扇形束重建問題。

FDK算法先對(duì)投影數(shù)據(jù)進(jìn)行預(yù)加權(quán)，適當(dāng)修正體素到坐標(biāo)系原點(diǎn)的距離和角度差；再對(duì)不同角度的投影數(shù)據(jù)進(jìn)行水平方向的一維濾波；最后沿射線方向進(jìn)行三維反投影。該過程可表示為

中間重建圖像與金屬圖像的融合過程為

其中，metal()為金屬圖像；metal為金屬區(qū)域；cor()為最終的校正圖像。

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

2.1 圖像質(zhì)量評(píng)價(jià)指標(biāo)

采用均方根誤差(root mean squared error，RMSE)、峰值信噪比(peak signal to noise ratio，PSNR)對(duì)校正后的重建圖像質(zhì)量進(jìn)行定量評(píng)價(jià)[22]。

RMSE定義為

PSNR定義為

PSNR可用來判斷圖像噪聲抑制能力，值越大表示抑制效果越好。PSNR是一種普遍采用的評(píng)估圖像質(zhì)量的客觀度量，有時(shí)無法與人的主觀感受完全匹配。

2.2 Shepp-Logan頭部模型數(shù)據(jù)校正實(shí)驗(yàn)

為了驗(yàn)證本文算法的有效性，在原三維Shepp-Logan頭部模型[23]中加入2個(gè)金屬球，進(jìn)行金屬偽影校正實(shí)驗(yàn)，模型參數(shù)見表1。

表1 三維Shepp-Logan頭部模型參數(shù)[23]

表1中，為橢球的附加強(qiáng)度值，為橢球的中心坐標(biāo)，，，分別為橢球在，，軸方向上的半軸長(zhǎng)系數(shù)(長(zhǎng)度為130 mm)，，，分別表示橢球繞，，軸旋轉(zhuǎn)的角度；序號(hào)1~10為原Shepp-Logan模型參數(shù)，序號(hào)11~12為金屬球參數(shù)。Shepp-Logan頭部模型中的多個(gè)橢球組成不同類型的組織，包括空氣、脂肪、骨、軟組織和水。三維Shepp-Logan頭部模型尺寸為256×256×256、體素尺寸為0.5859×0.5859×0.5859 mm3(各向同性)，圖2給出模型的向第128層切片、向第115層切片和向第185層切片。

采用CBCT掃描方式，對(duì)三維Shepp-Logan頭部模型進(jìn)行正向投影。射線源到平板探測(cè)器中心距離為1 000 mm，到旋轉(zhuǎn)中心的距離為550 mm；射線源的管電壓為140 kV，射線能譜采用BOONE和SEIBERT[24]提出的TASMIP算法得到；平板探測(cè)器的像素矩陣為560×560，平板探測(cè)器的像元尺寸為0.5×0.5 mm2；掃描范圍為0~360°，角度步長(zhǎng)為1°，投影角共360個(gè)。

在40~140 kV的射線能譜上，取能量步長(zhǎng)為10 kV，在美國(guó)標(biāo)準(zhǔn)與技術(shù)研究院(National Institute of Standards and Technology，NIST)官網(wǎng)上[5]，分別查找空氣、脂肪、骨、軟組織和水、鈦合金對(duì)應(yīng)不同射線能量的吸收系數(shù)。利用射線驅(qū)動(dòng)模型對(duì)各能量段的吸收系數(shù)進(jìn)行正向投影[25]，獲得單能投影數(shù)據(jù)，再經(jīng)加權(quán)融合，獲得最終的投影數(shù)據(jù)[26]，即

其中，w為第個(gè)能量段的權(quán)重系數(shù)；為第個(gè)角度上第個(gè)能量段的投影數(shù)據(jù)；為能量段總數(shù)，值為10；p為第個(gè)角度的投影數(shù)據(jù)，圖3給出第90個(gè)角度(90°角)對(duì)應(yīng)的投影數(shù)據(jù)。

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

圖2 三維Shepp-Logan頭部模型。圖像尺寸為256×256，灰度窗口為[0, 1]

圖3 Shepp-Logan模型的90°角投影數(shù)據(jù)。圖像尺寸為560×560，灰度窗口為[0.0001, 3.3128]

采用本文PIB-MAR算法進(jìn)行校正，圖5為該算法的各階段結(jié)果。圖5(a)為BF處理后的圖像，原圖像中的噪聲和輕微偽影得到平滑，邊緣信息得到保留；圖5(b)為金屬閾值分割圖像，金屬得到準(zhǔn)確分割，邊緣清晰；圖5(c)為先驗(yàn)圖像，各組織得到準(zhǔn)確聚類；圖5(d)為投影修復(fù)圖像，金屬投影區(qū)域連續(xù)性較好(對(duì)比圖3)；圖5(e)為中間重建圖像，原金屬區(qū)域連續(xù)性較好；圖5(f)為最終的校正圖像，金屬偽影得到有效校正(對(duì)比圖4(a))。

分別采用常規(guī)的LIB-MAR和IIB-MAR算法對(duì)原重建圖像進(jìn)行校正，并與本文的PIB-MAR算法進(jìn)行對(duì)比，結(jié)果如圖6所示。

圖5 基于先驗(yàn)圖像的Shepp-Logan模型的金屬偽影校正各階段結(jié)果

圖6 Shepp-Logan模型的金屬偽影校正結(jié)果對(duì)比。第1~3列分別為Z向第128層切片、Y向第115層切片和X向第185層切片；第1~3行分別為L(zhǎng)IB-MAR，IIB-MAR和PIB-MAR算法校正結(jié)果。圖像尺寸為256×256，灰度窗口為[0, 0.07]

從圖6中可以看出，相比原重建圖像(圖4)，LIB-MAR算法對(duì)應(yīng)的校正圖像偽影有所減少，兩金屬間的暗帶也被消除，但引起了次級(jí)偽影，部分組織結(jié)構(gòu)被破壞；與LIB-MAR算法相比，IIB-MAR算法在組織結(jié)構(gòu)保護(hù)上有更好的效果，但引起了更嚴(yán)重的次級(jí)偽影；與前2種MAR算法相比，PIB-MAR算法對(duì)應(yīng)的校正圖像偽影較少且組織結(jié)構(gòu)較完整，與原始圖像(圖2)最接近，獲得了最優(yōu)的校正效果。

為了直觀地對(duì)比3種算法的校正效果，繪制向第128層切片第116行的灰度曲線，如圖7所示。其中，參考灰度曲線來自原Shepp-Logan圖像 (圖2)。

從圖7(b)~(d)可以看出，PIB-MAR算法對(duì)應(yīng)的紅色虛線與參考圖像對(duì)應(yīng)的綠線最貼近，這表明該算法對(duì)應(yīng)的校正圖像質(zhì)量最優(yōu)。

為了定量地比較3種MAR算法的優(yōu)劣，在偽影較為嚴(yán)重的向第128層切片上選取3個(gè)尺寸為25×30的感興趣區(qū)域(region of interest，RoI)，如圖6所示。分別在單個(gè)及所有RoI上計(jì)算不同算法對(duì)應(yīng)的RMSE值和PSNR值，結(jié)果見表2。其中，參考圖像為原Shepp-Logan圖像(圖2)。

圖7 Shepp-Logan模型的金屬偽影校正結(jié)果的Z向第128層切片第116行的灰度曲線 ((a)為完整的灰度曲線; (b)~(d)分別是(a)中黑色虛框1~3的放大圖)

表2 單個(gè)及所有RoI區(qū)域上，不同MAR算法對(duì)應(yīng)的RMSE值和PSNR值

表2中，在同一RoI上，PIB-MAR算法對(duì)應(yīng)的RMSE值最??；在所有RoI上，PIB-MAR算法對(duì)應(yīng)的RMSE為0.002 1，比LIB-MAR和IIB-MAR算法分別減少了12.50%和41.67%。表明PIB-MAR算法對(duì)應(yīng)的校正圖像與原始圖像的偏差最小。在同一RoI中，PIB-MAR算法對(duì)應(yīng)的PSNR值最大；在所有RoI上，PIB-MAR算法對(duì)應(yīng)的PSNR為53.439 1，比LIB-MAR和IIB-MAR算法分別提高了1.65%和7.93%。表明該算法對(duì)校正圖像中的噪聲具有最強(qiáng)的抑制能力。

2.3 臨床頭部數(shù)據(jù)校正實(shí)驗(yàn)

為了進(jìn)一步驗(yàn)證PIB-MAR算法的性能，采用臨床頭部數(shù)據(jù)進(jìn)行MAR實(shí)驗(yàn)。臨床數(shù)據(jù)經(jīng)降采樣、歸一化、金屬植入處理，圖像尺寸為256×256×138 mm3，體素尺寸為0.5 mm，層厚為1 mm，中心切片如圖8(a)所示。采用上述CBCT掃描方式，對(duì)臨床數(shù)據(jù)進(jìn)行正向投影、FDK重建，得到原重建圖像，如圖8(b)所示，其中可見明顯的金屬偽影。分別采用LIB-MAR，IIB-MAR和PIB-MAR算法對(duì)原重建圖像進(jìn)行MAR，結(jié)果如圖8(d)~(f)所示。

由圖8可知，LIB-MAR算法去除了金屬周圍的大部分偽影，也基本去除了右下角的條紋狀偽影，但校正圖像中仍存在較多的次級(jí)偽影；IIB-MAR算法雖然去除了暗帶狀金屬偽影，但圖像中金屬鄰域連續(xù)性差且存在嚴(yán)重的次級(jí)偽影。與2種MAR算法相比，PIB-MAR算法去除了絕大部分金屬偽影，校正圖像中次級(jí)偽影輕微，取得了最佳的金屬偽影校正效果。

為了定量比較3種算法的優(yōu)劣，選取3個(gè)尺寸為30×30的RoI進(jìn)行RMSE值、PSNR值計(jì)算，結(jié)果見表3。

表3中，在同一RoI上，PIB-MAR算法對(duì)應(yīng)的RMSE值最小；在所有RoI上，其對(duì)應(yīng)的RMSE為0.002 4，與LIB-MAR和IIB-MAR算法相比，分別減少了31.42%和69.23%。這表明PIB-MAR 算法對(duì)應(yīng)的校正圖像與原始圖像的偏差最小。在同一RoI上，PIB-MAR算法對(duì)應(yīng)的PSNR值最大；在所有RoI上，其對(duì)應(yīng)的PSNR為57.122 4，與LIB-MAR和IIB-MAR算法相比，分別提高了18.15%和42.44%。這表明PIB-MAR算法具有最強(qiáng)的金屬偽影抑制能力。

圖8 臨床數(shù)據(jù)的金屬偽影校正結(jié)果。圖像尺寸為256×256，灰度窗口為[0, 0.05]

表3 Shepp-Logan體模的單個(gè)及所有RoI上，不同MAR算法對(duì)應(yīng)的RMSE和PSNR值

3 結(jié) 論

本文提出了一種基于先驗(yàn)圖像的CBCT金屬偽影校正算法，該算法首先對(duì)含金屬偽影的重建圖像雙邊濾波、金屬閾值分割、組織聚類等預(yù)處理，可較好地去除重建圖像中的噪聲并保留邊緣信息，從而獲得金屬圖像和不含金屬信息的先驗(yàn)圖像；再對(duì)二者正向投影，獲得金屬投影區(qū)域和先驗(yàn)投影數(shù)據(jù)；而后利用先驗(yàn)投影數(shù)據(jù)及金屬邊界鄰域的投影數(shù)據(jù)對(duì)金屬投影區(qū)域插值，獲得修復(fù)的投影數(shù)據(jù)，目的是抑制次級(jí)偽影產(chǎn)生；最后利用FDK算法對(duì)修復(fù)的投影數(shù)據(jù)重建，并與金屬圖像融合，獲得最終的校正圖像。利用三維Shepp-Logan頭部模型數(shù)據(jù)和臨床頭部CT數(shù)據(jù)進(jìn)行金屬偽影校正實(shí)驗(yàn)，結(jié)果表明，該算法在RMSE，PSNR 2項(xiàng)指標(biāo)上，均優(yōu)于常用的LIB-MAR和IIB-MAR算法，獲得較好的金屬偽影校正效果。

[1] GOVILA S, GUNDAPPA M. Cone beam computed tomography - an overview[J]. Journal of Conservative Dentistry, 2007, 10(2): 53-58.

[2] VENKATESH E, ELLURU S V. Cone beam computed tomography: basics and applications in dentistry[J]. Journal of Istanbul University Faculty of Dentistry, 2017, 51(3): 102-121.

[3] BORNSTEIN M, SCARFE W, VAUGHN V, et al. Cone beam computed tomography in implant dentistry: a systematic review focusing on guidelines, indications, and radiation dose risks[J]. The International Journal of Oral & Maxillofacial Implants, 2014, 29(Supplement): 55-77.

[4] 肖文, 曾理. CT圖像的金屬偽影校正方法綜述[J]. 中國(guó)體視學(xué)與圖像分析, 2019, 24(1): 29-36. XIAO W, ZENG L. The review of metal artifact reduction for CT images[J]. Chinese Journal of Stereology and Image Analysis, 2019, 24(1): 29-36 (in Chinese).

[5] 李銘. CT金屬偽影校正研究[D]. 北京: 中國(guó)科學(xué)院大學(xué), 2015. LI M. The research on metal artifacts reduction in computed tomography[D]. Beijing: University of Chinese Academy of Sciences, 2015 (in Chinese).

[6] WANG G, FREI T, VANNIER M W. Fast iterative algorithm for metal artifact reduction in X-ray CT[J]. Academic Radiology, 2000, 7(8): 607-614.

[7] ZHU Y C, LIU Y C, ZHANG Q, et al. A fast iteration approach to undersampled cone-beam CT reconstruction[J]. Journal of X-Ray Science and Technology, 2019, 27(1): 111-129.

[8] MEILINGER M, SCHUTZ O, SCHMIDGUNST C, et al. Alignment correction during metal artifact reduction for CBCT using mutual information and edge filtering[C]//2009 Proceedings of 6th International Symposium on Image and Signal Processing and Analysis. New York: IEEE Press, 2009: 135-140.

[9] ZHANG Y B, ZHANG L F, ZHU X R, et al. Reducing metal artifacts in cone-beam CT images by preprocessing projection data[J]. International Journal of Radiation Oncology·Biology·Physics, 2007, 67(3): 924-932.

[10] BEISTER M, PRELL D, KYRIAKOU Y, et al. GPU-accelerated metal artifact reduction (MAR) in FD-CT[C]//Medical Imaging 2010: Physics of Medical Imaging. San Diego: SPIE Press, 2010: 76223D.

[11] PRELL D, KYRIAKOU Y, STRUFFERT T, et al. Metal artifact reduction for clipping and coiling in interventional C-arm CT[J]. American Journal of Neuroradiology, 2010, 31(4): 634-639.

[12] BORNEMANN F, M?RZ T. Fast image inpainting based on coherence transport[J]. Journal of Mathematical Imaging and Vision, 2007, 28(3): 259-278.

[13] 張意, 蒲亦非, 張衛(wèi)華, 等. 基于一致性擴(kuò)散的非迭代快速CT金屬偽影去除算法[J]. 四川大學(xué)學(xué)報(bào): 工程科學(xué)版, 2013, 45(2): 81-87, 93. ZHANG Y, PU Y F, ZHANG W H, et al. A novel non-iterative metal artifact reduction algorithm using coherence transport in computed tomography[J]. Journal of Sichuan University: Engineering Science Edition, 2013, 45(2): 81-87, 93 (in Chinese).

[14] GUNTURK B. Fast bilateral filter with arbitrary range and domain kernels[C]//2010 IEEE International Conference on Image Processing. New York: IEEE Press, 2010: 3289-3292.

[15] OZA S, JOSHI K R. CUDA based fast bilateral filter for medical imaging[C]//2018 5th International Conference on Signal Processing and Integrated Networks (SPIN). New York: IEEE Press, 2018: 930-935.

[16] PAPARI G, IDOWU N, VARSLOT T. Fast bilateral filtering for denoising large 3D images[J]. IEEE Transactions on Image Processing, 2017, 26(1): 251-261.

[17] 陳金位, 吳冰. 二維直方圖重建和降維的Otsu閾值分割算法[J]. 圖學(xué)學(xué)報(bào), 2015, 36(4): 570-575. CHEN J W, WU B. A otsu threshold segmentation method based on rebuilding and dimension reduction of the two-dimensional histogram[J]. Journal of Graphics, 2015, 36(4): 570-575 (in Chinese).

[18] 陳培興, 王明泉, 李世虎, 等. 基于形態(tài)學(xué)和Otsu的固體火箭發(fā)動(dòng)機(jī)CT缺陷三維分割[J]. 圖學(xué)學(xué)報(bào), 2015, 36(4): 581-586. CHEN P X, WANG M Q, LI S H, et al. 3D CT defects segmentation of solid rocket motor based on the morphology and otsu[J]. Journal of Graphics, 2015, 36(4): 581-586 (in Chinese).

[19] BAL M, SPIES L. Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[J]. Medical Physics, 2006, 33(8): 2852-2859.

[20] SHI B Q, LIANG J, LIU Q. Adaptive simplification of point cloud using k-means clustering[J]. Computer-Aided Design, 2011, 43(8): 910-922.

[21] RODET T, NOO F, DEFRISE M. The cone-beam algorithm of Feldkamp, Davis, and Kress preserves oblique line integrals[J]. Medical Physics, 2004, 31(7): 1972-1975.

[22] LIU Y C, ZHANG C S, LI W R, et al. An adaptive multiscale anisotropic diffusion regularized image reconstruction method for digital breast tomosynthesis[J]. Australasian Physical & Engineering Sciences in Medicine, 2018, 41(4): 993-1008.

[23] 劉澤, 孫豐榮, 李艷玲, 等. 基于3D Shepp-Logan頭部模型的三維醫(yī)學(xué)圖像重建仿真[J]. 生物醫(yī)學(xué)工程學(xué)雜志, 2006(5): 938-943. LIU Z, SUN F R, LI Y L, et al. The three-dimension medical image reconstruction simulation on 3D Shepp-Logan head phantom[J]. Journal of Biomedical Engineering, 2006, 23(5): 938-943 (in Chinese).

[24] BOONE J M, SEIBERT J A. An accurate method for computer-generating tungsten anode X-ray spectra from 30 to 140 kV[J]. Medical Physics, 1997, 24(11): 1661-1670.

[25] SIDDON R L. Fast calculation of the exact radiological path for a three-dimensional CT array[J]. Medical Physics, 1985, 12(2): 252-255.

[26] 陳平, 潘晉孝, 劉賓. 連續(xù)能譜X-CT投影仿真算法[J]. 無損檢測(cè), 2009, 31(2): 102-104. CHEN P, PAN J X, LIU B. Simulation arithmetic of X-CT projection based on consecutive spectrum[J]. Nondestructive Testing, 2009, 31(2): 102-104 (in Chinese).

A prior-image-based metal artifact reduction method for cone beam CT

LIU Yang-chuan1, GAO Peng1,2, ZHU Ye-chen1, GAO Xin1

(1. Medical Imaging Department, Suzhou Institute of Biomedical Engineering and Technology, Chinese Academy of Sciences, Suzhou Jiangsu 215163, China; 2. School of Electronic and Optical Engineering, Nanjing University of Science and Technology, Nanjing Jiangsu 210094, China)

To effectively suppress the artifacts caused by metal implants in the reconstruction process of cone beam CT (CBCT) image, a prior-image-based metal artifact reduction method was proposed. Firstly, the reconstructed image with metal artifacts was preprocessed by bilateral filtering, metal threshold segmentation and tissue clustering to produce the metal image and the prior image without metal information. Secondly, the metal image and prior image were respectively forward-projected to produce the metal projection region and prior projection data. Then, the metal projection region was interpolated by the prior projection data and the metal neighborhood projection data to produce the restored projection data. Finally, the CT image was reconstructed by the FDK algorithm and was fused with the metal image to produce the final corrected image. To verify the performance of the proposed algorithm, the metal artifact reduction experiments were carried out on the 3D Shepp-Logan head phantom and clinical head CT data. The experimental results show that compared with the commonly used linear-interpolation-based method and image-inpainting-based method, the corrected image of the proposed method can keep the root-mean-square error to the minimum and the peak signal-to-noise ratio to the maximum. This indicates that the proposed method can effectively suppress metal artifacts while preserving image edge information.

cone beam CT; metal artifact reduction; bilateral filtering; prior image; interpolation

TP 751.1

10.11996/JG.j.2095-302X.2020040529

A

2095-302X(2020)04-0529-10

2020-02-11；

2020-04-22

22 April，2020

11 February，2020；

國(guó)家自然科學(xué)基金項(xiàng)目(61801475)；中國(guó)博士后科學(xué)基金項(xiàng)目(2018M642320)；江蘇省博士后科研資助項(xiàng)目(2018K180C)；中科院蘇州醫(yī)工所自主部署項(xiàng)目(Y95K091K05)；天津市科技計(jì)劃項(xiàng)目(19YDYGHZ00030)

National Natural Science Foundation of China (61801475); China Postdoctoral Science Foundation (2018M642320); Jiangsu Planned Projects for Postdoctoral Research Funds (2018K180C); Suzhou Institute of Biomedical Engineering and Technology (CAS) Planned Projects (Y95K091K05); Science and Technology Plan Projects of Tianjing (19YDYGHZ00030)

劉仰川(1987-)，男，山東濟(jì)寧人，副研究員，博士。主要研究方向?yàn)閄射線斷層成像技術(shù)。E-mail：liuyangchuan2006@163.com

LIU Yang-chuan (1987-), male, associate professor, Ph.D. His main research interests cover X-ray tomography. E-mail：liuyangchuan2006@163.com

高 欣(1975-)，男，吉林吉林人，研究員，博士。主要研究方向?yàn)榻槿朐\療。E-mail：xingaosam@yahoo.com

GAO Xin (1975-), male, professor, Ph.D. His main research interests cover interventional diagnosis and treatment. E-mail：xingaosam@yahoo.com