李 銘，盧彥飛，袁 剛，吳中毅，張 濤

（1.中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇蘇州215163；

2.中國科學(xué)院長春光學(xué)精密機械與物理研究所，吉林長春130033；

3.中國科學(xué)院大學(xué)，北京100049）

應(yīng)用先驗插值校正CT金屬偽影

李 銘1，2，3?，盧彥飛1，2，3，袁 剛1，吳中毅1，張 濤2

（1.中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇蘇州215163；

2.中國科學(xué)院長春光學(xué)精密機械與物理研究所，吉林長春130033；

3.中國科學(xué)院大學(xué)，北京100049）

針對CT系統(tǒng)在實際應(yīng)用中出現(xiàn)的金屬偽影問題，提出一種基于先驗插值的金屬偽影校正算法.文中通過預(yù)濾波、骨骼分割和軟組織恢復(fù)步驟計算先驗圖像，并利用先驗圖像的正向投影對原始投影中的金屬投影區(qū)進行插值校正.應(yīng)用該算法對數(shù)值仿真圖像和臨床CT圖像分別進行了校正重建實驗.數(shù)值仿真實驗表明，用提出算法校正的結(jié)果比線性插值金屬偽影校正算法、歸一化金屬偽影校正算法校正的結(jié)果更接近理想體模.臨床數(shù)據(jù)實驗表明:該算法的重建結(jié)果有效抑制了金屬偽影，清晰重建出金屬邊緣細節(jié)，極大地提高了重建圖像的質(zhì)量.

金屬偽影；先驗圖像；預(yù)濾波；軟組織恢復(fù)

1　引 言

計算機斷層成像（Computed Tomography，CT）是通過無損方式獲取人體內(nèi)部解剖信息的一種醫(yī)學(xué)影像技術(shù).CT掃描中，植入病人體內(nèi)的金屬物，如:股骨內(nèi)固定、金屬假牙、頸椎內(nèi)固定等可以在重建結(jié)果中產(chǎn)生偽影.這些偽影降低了斷層圖像中金屬周圍的清晰度，直接影響醫(yī)生對斷層結(jié)構(gòu)的判斷.

金屬偽影［1-2］主要表現(xiàn)為由大投影數(shù)據(jù)引起的從金屬區(qū)域發(fā)出的條狀偽影，由X射線能譜硬化引起的杯狀偽影和多個金屬之間的暗帶區(qū)域.金屬偽影使重建圖像的質(zhì)量嚴重下降，因此需要有效的方法從圖像域或投影域?qū)ζ溥M行校正.金屬偽影校正算法較多，可以大致分為兩類:投影插值法［3-8］和迭代法［9-12］.迭代法能夠有效處理復(fù)雜的偽影問題，但由于計算量過大而限制其使用.近年來，基于投影空間插值算法形成了金屬分割、正向投影和投影插值的較為統(tǒng)一的校正框架. Gu等［5］采用線性插值對原始投影空間進行插值處理.線性插值算法（LI-MAR）能有效地消除金屬周圍的暗影區(qū)域，卻容易模糊金屬邊緣的組織結(jié)構(gòu)和引入新的正切于金屬邊緣的次級偽影. Bal等［6］引入均值聚類先驗圖像來校正投影空間中的金屬投影區(qū).該算法能有效地保護金屬邊緣的骨骼結(jié)構(gòu)，但對于次級偽影的抑制效果不明顯. Lell等［8］針對Bal算法的不足，應(yīng)用均值聚類先驗圖像做標(biāo)準化插值，提出了歸一化金屬偽影校正（NMAR）算法.NMAR算法能更有效地抑制次級偽影和保護金屬周圍的高密度組織.

本文研究仍采用上述的算法框架，針對均值聚類先驗圖像具有軟組織信息不完整的缺陷，提出了預(yù)濾波、骨骼分割和軟組織恢復(fù)處理的新的先驗圖像計算方法，并在此基礎(chǔ)上提出一種基于先驗插值的金屬偽影校正（PI-MAR）算法.采用體模仿真圖像和真實CT圖像對新算法進行了實驗，并對其與LI-MAR算法和NMAR算法進行了比較.結(jié)果表明，所提出的算法能更有效地去除CT圖像中的金屬偽影和保護斷層圖像中金屬周圍的組織結(jié)構(gòu).

2　先驗圖像計算

N為局部濾波窗口半徑，w（i，j）為依賴于中心像素i和鄰域像素j相似程度的權(quán)重.定義中心像素i和鄰域像素j間的距離越大，則相似程度越低；若兩者之間的CT值超過給定閾值，則定義它們間的權(quán)重為0，即鄰域像素j不參與中心像素i的計算.

式中:Z（i）為歸一化因子，通過式（3）可得:

2.1預(yù)濾波處理

為了抑制原始圖像中的條紋狀偽影和保護組織結(jié)構(gòu)，本文提出了帶閾值約束的局部均值濾波算法.該算法的基本思想是:當(dāng)前位置的像素由局部窗口內(nèi)所有與它相似的像素加權(quán)平均得到.區(qū)域生長圖像分割算法能將具有相同特征的連通域分割出來，同時能提供良好的邊界信息，而金屬物區(qū)域通常為連通區(qū)域，金屬區(qū)域和普通組織又具有明顯的邊界，因此本文采用基于區(qū)域生長的圖像分割算法提取原始圖像中的金屬物.具體執(zhí)行預(yù)濾波處理前，首先把金屬物對應(yīng)位置的像素其中:T為自定義閾值，用以識別組織像素和偽影像素；h為濾波參數(shù)，用以控制濾波強度.實際計算過程當(dāng)中，關(guān)鍵濾波參數(shù)T的設(shè)置至關(guān)重要.由于條狀偽影和軟組織之間的像素值差通常會大于80 HU，為了更有效地抑制條狀偽影和保護組織結(jié)構(gòu)，文中選擇的閾值T小于80 HU.具體實驗中:N=5，T=60 HU，h=10能取得較好的偽影抑制和組織邊緣保護效果.

2.2骨骼分割

依據(jù)預(yù)濾波處理后圖像的聚類特性，本文將Mean shift算法應(yīng)用于濾波后的圖像，提取骨骼結(jié)構(gòu).對比閾值分割算法，Mean shift算法能更精確地提取骨骼結(jié)構(gòu)信息.當(dāng)Mean shift算法［13-14］應(yīng)用于圖像分割時，統(tǒng)一考慮圖像的空間信息和灰度信息，組成一個p+2維的向量x=（xs，xr），其中xs表示網(wǎng)格點的坐標(biāo)，xr表示該網(wǎng)格點上的p維向量特征.當(dāng)圖像為灰度圖像時，p=1；彩色RGB圖像時，p=3.分別用xi和zi，i=1，…，n表示原始和分割后圖像.定義向量

式中:w（xi）為采樣點x的權(quán)重，離x近的采樣點xi有較大的權(quán)重，即離x越近的采樣點，對估計x周圍的統(tǒng)計特性越有效，反之亦然.

根據(jù)參考文獻［14］，定義核函數(shù)G（x）具有如下形式:

其中:C為歸一化常數(shù)，hs，hr為核帶寬.hs為空間帶寬，其值越大，分割需要時間越長；hr為像素點的灰度帶寬，其值越大，越多的圖像細節(jié)就會被忽略.為了保證一定的計算速度和分割分辨率，hs，hr均不宜選取特別大的值.

用Mean shift算法進行圖像分割的具體步驟如下（對每一像素點）:

（1）初始化l=1，并且使yi，l=xi

（2）依據(jù)公式（4）計算yi，l+1=MS（yi，l）

（3）重復(fù)步驟（2）直到y(tǒng)i，l收斂到最終結(jié)果yi，L，則輸出結(jié)果zi=（xsi，yri，L）

（4）將收斂至同一起始點的歸為一類，合并像素點過少的類，得到最終結(jié)果.

結(jié)合Mean shift算法處理結(jié)果，采用簡單閾值即可提取濾波后CT圖像的骨骼結(jié)構(gòu).

2.3軟組織恢復(fù)

植入病人體內(nèi)的高密度物質(zhì)，可以在重建圖像的金屬周圍產(chǎn)生亮暗區(qū)域.這些亮暗偽影的存在，嚴重影響了對斷層結(jié)構(gòu)的判斷，使圖像質(zhì)量下降，難以達到無損檢測、醫(yī)療診斷等目的.因此，需要恢復(fù)靠近金屬周圍的軟組織結(jié)構(gòu).為了恢復(fù)靠近金屬鄰域的軟組織信息，首先應(yīng)識別金屬亮暗偽影的具體位置.本文通過在偽影區(qū)域設(shè)置種子點，使用區(qū)域生長方法識別靠近金屬位置的偽影像素.為了使恢復(fù)的像素值能夠反映軟組織結(jié)構(gòu)的浮動，并且與周圍的像素具有相對平滑的過渡，文中構(gòu)造如下的軟組織賦值函數(shù):

式中:D表示賦值點與距離最近的金屬像素的歐氏距離，Imin表示偽影區(qū)域邊界處的最小值，c是曲線參數(shù)，A表示軟組織CT值的浮動范圍.通常軟組織的浮動范圍會小于150 HU，為了保證賦值的像素與周圍的像素具有相對平滑的過渡，下文的具體實驗中選擇A=150 HU.使用上述的構(gòu)造函數(shù)恢復(fù)偽影區(qū)域的像素值后，將2.2節(jié)提取的骨骼結(jié)構(gòu)植入，得到最終的先驗圖像.

3　金屬投影區(qū)插值

最后，使用濾波反投影法對校正后的投影數(shù)據(jù)進行重建，并植入金屬物，得到最終的校正重建圖像.本文提出的金屬偽影校正流程如圖（1）所示.

圖1　金屬偽影校正流程Fig.1　Flowchart of the proposed MAR method

4　實驗與結(jié)果

為了證明本文算法在抑制CT金屬偽影和保護金屬周圍組織方面的有效性，分別采用體模仿真圖像和臨床CT圖像對本文提出的算法進行實驗驗證.實驗采用CT圖像的分辨率為512× 512，投影數(shù)據(jù)的分辨率為1 024×720.實驗平臺:Intel（R）Core（TM）2 Quad CPU Q8200 2.66 GHz/內(nèi)存4 GB，實驗程序采用Matlab編寫實現(xiàn).

4.1仿真實驗

本部分實驗選擇體模仿真圖像分別對LIMAR、NMAR和PI-MAR算法進行實驗驗證.具體數(shù)值仿真實驗中:CT射線源的能譜采用Boone等［16］提出的物理模型模擬得到，圖2所示為文中用于仿真實驗的射線能譜圖；本文采用Poisson分布模擬探測器接受X射線光子這一隨機過程，具體數(shù)學(xué)模型［17］如下所示:

式中:gi為探測器i實際接受的光子數(shù)，為k能量下入射光子數(shù)，Li為X射線的路徑，μ表示k能量下組織的衰減系數(shù)，Sik用于統(tǒng)計k能量下的散射光子數(shù)，K表示對X射線源能譜的實際采樣數(shù).通過對公式（9）計算的結(jié)果，做對數(shù)化處理后，生成X射束硬化的投影數(shù)據(jù)，緊接著對投影數(shù)據(jù)做反投影重建得到仿真CT圖像.理想體模圖像由不含金屬的數(shù)值模型仿真，再植入金屬物得到；原始圖像由含金屬的數(shù)值模型仿真得到.仿真實驗采用理想體模圖像作為3種算法校正的結(jié)果對比的金標(biāo)準.具體仿真實驗結(jié)果如圖3所示，由圖3可以看到，仿真原始圖像的金屬周圍出現(xiàn)嚴重的條狀偽影和射束硬化偽影，總體圖像質(zhì)量很差.采用LIMAR算法對其校正后，圖像質(zhì)量有了部分改善，但校正后圖像中仍存在明顯的帶狀偽影. NMAR算法校正的圖像對比LI-MAR算法，能取得更好的次級偽影抑制效果，總體圖像質(zhì)量比原始圖像也有了明顯的提升，但金屬邊緣仍有很少的次級偽影.而PI-MAR算法重建結(jié)果與理想體模圖像非常接近，金屬邊緣細節(jié)清晰，次級偽影抑制效果最好.

圖2　模擬120 ke V下球管的X射線能譜Fig.2　X-ray spectrum of 120 keV

圖3　仿真數(shù)據(jù)實驗Fig.3　Synthetic dataset experiment

為了定量評價上述3種金屬偽影校正算法的具體表現(xiàn)，從圖3（a）所示圖像位置，選定5個窗口區(qū)域，并計算感興趣區(qū)域（Region of interest，ROI）的CT均值和標(biāo)準差（Standard Deviation， SD）.ROI內(nèi)的CT均值反映校正后圖像與理想數(shù)值模型的逼真度，SD反映校正后的圖像ROI區(qū)域的數(shù)據(jù)波動，用來反映校正后圖像中次級偽影的強度.具體計算結(jié)果如表1所示.

表1　所有ROI位置的CT均值和標(biāo)準差對比Tab.1　Comparison of the mean CT number and the SD in all the ROIs

從表1可以看到，所有ROI位置，PI-MAR算法的CT均值與理想體模圖像的CT均值更接近，表明用提出的PI-MAR算法校正的圖像與金標(biāo)準的逼真度最高.與原始圖像對應(yīng)ROI的標(biāo)準差相比，通過3種算法校正后，幾乎所有ROI位置的標(biāo)準差都有不同程度的降低，這表明3種算法均能有效抑制金屬偽影.對比3種算法所有ROI位置的標(biāo)準差，PI-MAR算法校正結(jié)果的標(biāo)準差最小，表明PI-MAR算法對次級偽影的抑制效果更好.

圖4　臨床數(shù)據(jù)實驗Fig.4　Clinical dataset experiment

4.2臨床數(shù)據(jù)實驗

為了進一步驗證PI-MAR算法的性能，采用真實CT圖像作為原始實驗數(shù)據(jù)，分別對上述3種算法進行實驗.圖4給出臨床數(shù)據(jù)實驗結(jié)果，窗口顯示范圍為-500～500 HU.由圖4中原始圖像所示，金屬偽影在臨床CT圖像中具體表現(xiàn)為從金屬周圍發(fā)出的射線狀條紋和靠近金屬邊緣的低亮度區(qū)域.由于金屬偽影的存在，金屬邊緣的組織結(jié)構(gòu)很模糊，影響醫(yī)生對斷層組織結(jié)構(gòu)的判斷.使用LI-MAR算法對其校正后，靠近金屬周圍的暗偽影得到了有效抑制，但靠近金屬的骨骼模糊嚴重，次級偽影也很明顯.與LI-MAR算法的校正圖像相比，NMAR算法校正的結(jié)果對金屬周圍骨骼保護和偽影抑制均能取得更好的效果，但接近金屬的位置仍有少量的偽影存在.對比三種算法的校正結(jié)果可以看出，PI-MAR算法校正的圖像質(zhì)量最優(yōu)，能夠重建出清晰的金屬邊緣細節(jié)，偽影也得到了最有效抑制.然而在臨床數(shù)據(jù)實驗中，由于本文算法是從含有金屬偽影的CT圖像出發(fā)，進行金屬偽影校正，因此，校正后CT圖像分辨率會有所下降，若能直接從CT掃描設(shè)備上獲取相應(yīng)的投影數(shù)據(jù)，該問題便可迎刃而解.

為了更好地評價LI-MAR、NMAR和PIMAR算法校正結(jié)果的圖像質(zhì)量.我們讓兩位放射科專家采用Verburg等［18］提出的圖像質(zhì)量主管評價標(biāo)準對校正結(jié)果做評定.圖像質(zhì)量主觀評價結(jié)果如表2所示，原始圖像的平均圖像質(zhì)量等級為5.0，這表明原始圖像的質(zhì)量最差.通過上述3種算法分別對其校正后，其平均圖像質(zhì)量等級分別為3.0、1.9、1.0.這說明對比原始圖像，3種算法均能夠不同程度地消除金屬偽影，從而提升圖像質(zhì)量.從表2可以看出，PI-MAR算法校正的結(jié)果平均圖像質(zhì)量等級最低，表明本文提出的PI-MAR算法校正結(jié)果最優(yōu)，具有更好的臨床應(yīng)用前景.

表2　圖像質(zhì)量等級Tab.2　Rank of image quality

5　結(jié) 論

臨床CT診斷中金屬偽影的存在直接降低了CT圖像的質(zhì)量和其臨床診斷的準確性.本文提出一種基于先驗插值的金屬偽影校正算法.該算法通過預(yù)濾波、骨骼提取和軟組織恢復(fù)步驟計算具有完整組織信息的先驗圖像，解決了均值聚類先驗圖像軟組織信息不完整的缺陷.通過應(yīng)用一種平滑插值技術(shù)校正金屬投影區(qū)，有效抑制了次級偽影的形成，提高了校正圖像質(zhì)量.仿真實驗結(jié)果表明，文中算法重建圖像的CT均值較LIMAR算法、NMAR算法重建結(jié)果更接近理想的金標(biāo)準，證明本文PI-MAR算法重建的圖像與理想體模圖像更接近.臨床CT圖像實驗結(jié)果表明，PI-MAR算法有效抑制了金屬偽影，清晰重建出金屬邊緣細節(jié)，極大地提高了校正后的圖像質(zhì)量.

［1］ Robertson D D，Yuan J，Wang G，et al.Total hip prosthesis metal-artifact suppression using iterative deblurring reconstruction［J］.Journal of Computer Assisted Tomography，1997，21（2）:293-298.

［2］ 李銘，張濤，鄭健，等.基于切線反投影的CT金屬位置和形狀標(biāo)定［J］.液晶與顯示，2013，28（2）:295-299.

Li M，Zhang T，Zheng J，et al.Determination of location and shape of metallic object in CT based on tangent backprojection［J］.Chinese Journal of Liquid Crystals and Displays，2013，28（2）:295-299.（in Chinese）

［3］ Kalender W A，Hebel R，Ebersberger J.Reduction of CT artifacts caused by metallic implants［J］.Radiology，1987，164（2）:576-577.

［4］ 李印生，陳陽，馬建華，等.基于自適應(yīng)方向性濾波和非局部均值修補的CT圖像金屬偽影消除［J］.中國生物醫(yī)學(xué)工程學(xué)報，2011，30（3）:377-381.

Li Y S，Chen Y，Ma J H，et al.Metal artifact reduction in CT based on adaptive steering filter and nonlocal sinogram inpainting［J］.Chinese Journal of Biomedical Engineering，2011，30（3）:377-381.（in Chinese）

［5］ Gu J W，Zhang L，Chen Z Q，et al.A method based on interpolation for metal artifacts reduction in CT images［J］.Journal of X-Ray Science and Technology，2006，14（1）:11-19.

［6］ Bal M，Spies L.Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering［J］.Medical Physics，2006，33（8）:2852-2859.

［7］ Meyer E，Raupach R，Lell M，et al.Frequency split metal artifact reduction（FSMAR）in computed tomography［J］.Medical Physics，2012，39（4）:1904-1916.

［8］ Lell M M，Meyer E，Kuefner M A，et al.Normalized metal artifact reduction in head and neck computed tomography［J］.Investigative Radiology，2012，47（7）:415-421.

［9］ Oehler M，Buzug T M.Theλ-MLEM algorithm:An iterative reconstruction technique for metal artifact reduction in CT images［J］.Advances in Medical Engineering，2007，114:42-47.

［10］ Boas F E，F(xiàn)leischmann D.Evaluation of two iterative techniques for reducing metal artifacts in computed tomography［J］.Radiology，2011，259（3）:894-902.

［11］ Zhang X M，Xing L.Sequentially reweighted TV minimization for CT metal artifact reduction［J］.Medical Physics，2013，40（7）:071907.

［12］ Zhang Y B，Yan H，Jia X，et al.A hybrid metal artifact reduction algorithm for x-ray CT［J］.Medical Physics，2013，40（4）:041910.

［13］ 張宇洋，劉滿華，韓韜.基于Mean Shift圖像分割和支持向量機判決的候梯人數(shù)視覺檢測系統(tǒng)［J］.光學(xué)精密工程，2013，21（4）:1079-1085.

Zhang Y Y，Liu M H，Han T.Elevator-waiting people counting system based on Mean Shift segmentation and SVM classification［J］.Optics and Precision Engineering，2013，21（4）:1079-1085.（in Chinese）

［14］ Yu HY，Zeng K，Bharkhada D K，et al.A segmentation-based method for metal artifact reduction［J］.Academic Radiology，2007，14（4）:495-504.

［15］ Jacobs F，Sundermann E，De Sutter B，et al.A fast algorithm to calculate the exact radiological path through a pixel or voxel space［J］.Journal of Computing and Information Technology，1998，6（1）:89-94.

［16］ Boone J M，Seibert J A.An accurate method for computer-generating tungsten anode x-ray spectra from 30 to 140 k V［J］.Medical Physics，1997，24（11）:1661-1670.

［17］ La Rivière P J，Bian J G，Vargas P A.Penalized-Likelihood sinogram restoration for computed tomography［J］. IEEE Transactions on Medical Imaging，2006，25（8）:1022-1036.

［18］ Verburg J M，Seco J.CT metal artifact reduction method correcting for beam hardening and missing projections［J］.Physics in Medicine and Biology，2012，57（9）:2803-2818.

Metal artifact reduction in computed tomography based on prior interpolation

Li Ming1，2，3?，Lu Yan-fei1，2，3，Yuan Gang1，Wu Zhong-yi1，Zhang Tao2

（1.Suzhou Institute of Biomedical Engineering and Technology，Chinese Academy of Sciences，Suzhou 215163，China；
2.Changchun Institute of Optics，F(xiàn)ine Mechanics and Physics，Chinese Academy of Sciences，Changchun 130033，China；
3.University of Chinese Academy of Sciences，Beijing 100049，China）

Presence of metal artifacts in computed tomography（CT）is a major reason of degradation of image quality.In this work，a novel metal artifact reduction（MAR）algorithm using prior interpolation was developed to suppress metal artifacts in CT images.The prior image was generated by prefiltering，bone extraction，and soft tissue restoring procedures.This prior image was then reprojected，and the estimation of surrogate data in the metal traces could be performed via the reprojections.Next，the corrected experiments were conducted on digital phantom and clinical data respectively.The simulation experiment shows that results stem from prior interpolation MAR（PI-MAR）have higher similarity of ideal phantom than those from linear-interpolation MAR（LI-MAR）and Nor-malized MAR（NMAR）algorithms.Clinical experimental results demonstrate that the new algorithm can successfully suppress metal artifacts and clearly restore anatomical structure close to metallic implants.Moreover，the quality of reconstructed images can be greatly improved by the proposed algorithm.

metal artifacts；prior image；pre-filtering；soft tissue restoring

TP391.4

A doi:10.3788/YJYXS20153006.1032

1007-2780（2015）06-1032-08

李銘（1984-），男，河北石家莊人，博士研究生，主要從事CT圖像重建和金屬偽影校正方面的研究.E-mail: lim@sibet.ac.cn

2015-01-10；

2015-02-06.

國家自然科學(xué)基金（No.61301042，No.61201117）；國家重大科學(xué)儀器設(shè)備開發(fā)專項（No.2011YQ040082）；中國科學(xué)院青年創(chuàng)新促進會（No.2014281）；蘇州市科技計劃項目（No.ZXY2013001）；江蘇省自然基金項目（No.BK2012189，No.BK20151232）

Supported by National Natural Science Foundation of China（No.61301042，No.61201117）；Special Funds for the Development of National Major Scientific Instruments and Equipment（No.2011YQ040082）；Youth Innovation Promotion Association CAS（No.2014281）；Science and Technology Program of Suzhou（No. ZYX2013001）；Natural Science Foundation of Jiangsu Province（No.BK2012189，No.BK20151232）.

?通信聯(lián)系人，E-mail:lim@sibet.ac.cn