基于自適應(yīng)方向性濾波和非局部均值修補的CT圖像金屬偽影消除

李印生 陳 陽 馬建華 羅立民 陳武凡

1(東南大學(xué)計算機科學(xué)與工程學(xué)院，南京 210096)

2(南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，廣州 510515)

引言

CT系統(tǒng)中，許多因素可以在重建結(jié)果中產(chǎn)生偽影，在這些偽影中，金屬偽影最為明顯。CT掃描中，金屬等高衰減系數(shù)的物質(zhì)導(dǎo)致到達探測器的光子計數(shù)十分稀少，讀取的投影值就異常大，使得對于斷層結(jié)構(gòu)的判斷帶來很大困難。這些金屬物體包括外科手術(shù)中的金屬夾子、金屬修補物、填充物和金屬假牙等。由于這些高衰減系數(shù)的物質(zhì)大多都是各種金屬，所以由這些物體引起的偽影統(tǒng)稱為“金屬偽影”。

金屬偽影消除算法較多，可以大致分為三類：投影插值法、迭代法和混合方法。投影插值法理論簡單，實現(xiàn)快速，但對于全局放射狀偽影效果欠佳。迭代法能夠有效處理復(fù)雜的偽影問題，但由于計算量過大而限制其使用。混合法結(jié)合兩者優(yōu)點，目前受到廣泛關(guān)注。

近年來，基于投影(弦圖)空間插值算法形成了前置濾波、分割和投影插值的較為統(tǒng)一的框架[2]。前置濾波可以在一定程度上降噪并對于條狀偽影進行平滑處理。分割則是十分關(guān)鍵的一步，目的在于將圖像中的金屬成分和偽影成分與其他組織進行區(qū)分，繼而進行提取。分別對原始圖像和金屬成分進行前向投影，得到相應(yīng)的弦空間數(shù)據(jù)。做差運算即可將金屬偽影成分從弦空間中去除。然而，將金屬弦圖從原始圖像正弦圖中減除后，會造成原有的肌肉、軟組織等投影的不連續(xù)現(xiàn)象，反投影后，會產(chǎn)生新的偽影[2]。因此，需要對于投影空間數(shù)據(jù)進行插值處理。例如，2006年Matthieu等提出的金屬偽影消除算法就是采用了上述流程，取得一定的效果[2]。但是，Matthieu的算法有兩點明顯的不足：一是k-means聚類分割需要人工指定聚類數(shù)和聚類中心，兩者的設(shè)置隨圖像的不同需要做出調(diào)整，不具一般性;同時，聚類數(shù)的設(shè)置會影響到算法運行的效率和效果[4]。二是 Radon空間補全采用線性插值算法，利用局部信息補全弦圖中減除的部分。顯然，Radon空間補全還可以通過基于全局修補的方法以達到更好的結(jié)果。

本研究同樣采用上述算法框架，針對 Matthieu算法的不足，提出了基于方向性自適應(yīng)前置濾波，配合均值漂移(Mean shift)和最大互信息熵(MMS)分割[4]，最后使用非局部均值圖像修補(nonlocal means inpainting)[5]的新的金屬偽影消除算法。采用體模仿真圖像和真實CT圖像進行實驗，結(jié)果表明，所提出的算法能有效地去除CT圖像中的金屬偽影。

1 方法

1.1 自適應(yīng)方向性高斯濾波器

采用方向性自適應(yīng)濾波器(adaptive steering filter，ASF)，該濾波器具有濾波方向，尺度隨圖像局部幾何結(jié)構(gòu)而改變，并且濾波強度人工可調(diào);能夠很好的完成CT金屬偽影去除流程中前置濾波器的功能[2]。對于選取的各項異性高斯核闡述如下：

設(shè)G0為配分常數(shù)，σ1、σ2為各項異性高斯核的尺度系數(shù)，見圖1。各項異性高斯卷積核表述為

圖1 高斯卷積核平面投影Fig.1 The Gauss convolution kernel

式中，d1=u cosθ+v sinθ d2=u cosθ-v sinθ

其后，由非線性各向異性結(jié)構(gòu)張量(nonlinear anisotropic structure tensor，NAST)對圖像局部結(jié)構(gòu)進行估計，繼而確定ASF的尺度系數(shù)σ1、σ2和旋轉(zhuǎn)角度θ。結(jié)構(gòu)張量(以下簡稱 ST)通過求解圖像局部的梯度信息，構(gòu)造張量矩陣，根據(jù)其特征值和特征向量的具體形式判斷圖像的局部幾何結(jié)構(gòu)[8]。

首先定義圖像f的張量積矩陣

式中，fx、fy分別表示圖像對 x、y方向的一階偏導(dǎo)數(shù)。使用濾波算子 G(·)對張量矩陣進行濾波處理，當(dāng)G(·)為線性運算時，稱為線性結(jié)構(gòu)張量，當(dāng)G(·)為非線性運算時，稱為非線性結(jié)構(gòu)張量。線性結(jié)構(gòu)張量一般由式(2)中的張量矩陣與高斯核函數(shù)卷積得來，不可避免地造成圖像平滑模糊[8]。故采用非線性各向異性結(jié)構(gòu)張量估計圖像局部結(jié)構(gòu)，定義為

ε為補償系數(shù)，為了保證根號下的值為正數(shù)。p為邊緣強度，取 p=1.5進行實驗。最后，將 ?tTi，j對于擴散時間t進行積分，即可求得經(jīng)非線性擴散濾波處理后的NAST矩陣。

NSAT 矩陣特征值[8]與矩陣特征值對應(yīng)的特征向量可以反映圖像的局部結(jié)構(gòu)，例如，對于圖像均質(zhì)區(qū)域，灰度值變化程度比較平緩，張量矩陣的特征值較小或為0。對于圖像邊緣、角點或是被條狀偽影污染的部分，特征值較大[8]。此處ASF的作用就是要檢測出條狀偽影輻射的方向，并提供合適的卷積強度平滑條狀偽影。

根據(jù)以上闡述，首先確定偽影強度。ASF濾波器示意見圖2。對于 ASF的長軸尺度，采用文獻[2]中的方法進行估計。設(shè) σ1是條狀偽影參考坐標(biāo)系中與偽影正交方向的ASF濾波尺度，同樣，對于ASF的濾波方向，采用文獻[8]中的方法進行估計。式(1)中的G0為配分常數(shù)，使得

即可確定G0取值。

圖2 ASF濾波器示意Fig.2 ASF filter

最后，采用 ASF對圖像進行前置濾波。假設(shè)I0(x，y)為待處理圖像中某處灰度值，IASF(x，y)是經(jīng)過ASF處理后的同一位置處的圖像灰度值，Gx，y(u，v)為 ASF濾波核函數(shù)，則前置濾波過程可表述為

取r=2作為卷積核掩模的半徑。

1.2 最大互信息熵結(jié)合均值漂移的圖像金屬及偽影成分的分割

金屬及其偽影成分的精確分割是影響CT金屬偽影去除結(jié)果的重要因素。為了有效的分割出圖像中包含的金屬和條狀偽影，采用均值漂移結(jié)合最大互信息熵方法分割圖像中的金屬成分及其形成的條狀偽影。

最大互信息熵是針對k-均值聚類算法和Meanshift圖像分割算法這類區(qū)域生成算法的不足之處而專門設(shè)計的[4]。它基于模擬退火(Simulated Annealing)過程，以最大互信息量(MI)作為優(yōu)化分割目標(biāo)，以互信息熵差(d MI)作為分類數(shù)判據(jù)的一種優(yōu)化閾值分割算法。當(dāng)其小于某一閾值時，認(rèn)為分割結(jié)束。以原始圖像和當(dāng)前分類數(shù)下的分割圖像的互信息量為優(yōu)化目標(biāo)。當(dāng)互信息熵差達到最小值，同時互信息量達到最大值時，認(rèn)為達到最優(yōu)分割[4]?？梢?，該分割方法是典型的多聚類數(shù)分割，可以將金屬和條狀偽影及其他組織成分有效的區(qū)分。對于體模圖像和臨床真實圖像，選取的最大聚類數(shù)分別取為5和16。算法原理和實現(xiàn)過程具體參見文獻[4]。

1.3 結(jié)合線性插值的非局部均值Radon空間修補

在經(jīng)過前置濾波和分割去除偽影后，如果不對分割后的弦空間數(shù)據(jù)進行有效的修復(fù)，則做差運算后，代替原來金屬成分的部分與其周圍組織在弦空間的融合處不連續(xù)，重建后會產(chǎn)生新的偽影[2]。為了有效去除偽影，并最大限度地保證周圍其他組織不被破壞，本研究采用非局部均值修復(fù)的方法進行Radon空間的補全。基于較大遍歷窗的非局部均值修補，能夠有效利用圖像全局信息。圖像修復(fù)模型的算法原理為

Nonlocal權(quán)重系數(shù)為

Ni、Nj為以 i、j為中心，Rsim為半徑的一個小鄰域，h為衰減強度。

加權(quán)濾波定義為

Z歸一化常數(shù)定義為

聯(lián)立式(11)～(14)，得到用于圖像恢復(fù)的 Nonlocal加權(quán)項，完成圖像修補。

為了清楚闡述修補過程，定義以待修補區(qū)域像素為中心的鄰域為目標(biāo)窗。當(dāng)前計算權(quán)重的像素鄰域為相似窗。非局部均值修補的具體步驟如下：

首先以分割出金屬以及偽影成分弦圖數(shù)據(jù)定義掩模，將待修補圖像相應(yīng)位置用掩模覆蓋;繼而提取掩模部分(待修補區(qū)域)的邊界，并設(shè)置目標(biāo)窗、相似窗半徑Rsim，沿邊界順時針并逐步趨向待修補區(qū)域中心進行遍歷搜索;對于算法本身設(shè)置權(quán)重閾值δ和合適的衰減常數(shù)h;若當(dāng)前求得的相似權(quán)重w(i，j)≥δ，則將相似窗數(shù)據(jù)復(fù)制給目標(biāo)窗，否則，采用線性插值方法直接由目標(biāo)窗鄰近數(shù)據(jù)進行局部修補。

做出這樣的改進，是出于實際情況考慮的。首先，作者通過實驗發(fā)現(xiàn)，在有些情況下，無論如何設(shè)置窗口半徑Rsim和衰減常數(shù)h都可能無法找到與目標(biāo)窗相似程度較高的相似窗，這樣的修補無疑會產(chǎn)生較大誤差。此時應(yīng)采用線性插值由目標(biāo)窗附近的數(shù)據(jù)直接進行插值修補。并且，由待修補區(qū)域邊緣順時針向內(nèi)的修補也符合常用的圖像修補原則。實驗證明，這種方法是可行的。最后，只需將修補處理后的弦空間數(shù)據(jù)進行FBP重建，可得到本算法處理后的結(jié)果。

算法步驟為

步驟1：采用ASF前置濾波器對含有金屬偽影的圖像進行全局濾波;

步驟2：分別采用均值漂移分割濾波后圖像的金屬成分，采用最大互信息熵分割濾波后圖像的偽影成分;

步驟3：分別對包含偽影成分的圖像和原始圖像進行前向投影;

步驟4：在原圖投影中，依據(jù)金屬和偽影成分投影制作掩模;

步驟5：采用非局部均值結(jié)合線性插值方法對弦空間進行修補;

步驟6：采用FBP方法得到修補后弦空間數(shù)據(jù)的空間域圖像;

步驟7：將步驟6得到的圖像和均值漂移分割出的金屬成分進行空間域插值，得到最后的校正圖像。

1.4 算法驗證實驗

實驗分別采用體模和真實臨床數(shù)據(jù)，對本算法進行驗證和比較。CT圖像維數(shù)采用512像素×512像素。重建圖像采用 FBP算法，其中濾波核采用“B40 f”。校正算法中采用平行束成像幾何。圖3(a)是一個包含金屬偽影的模擬CT圖像，它包含了一個金屬植入物和一個圓柱形水容器。標(biāo)準(zhǔn)體模圖像基于下面的 CT值設(shè)置：空氣，0;方形水槽，500;圓柱形水槽，3000;水，1000。圖 3(b)是病人肝臟金屬縫合物產(chǎn)生的條狀偽影。計算機硬件環(huán)境為Inter Core i52.68×4 GHz，6 GB RAM。為了縮短校正時間，本研究采用混合Matlab/C++編程環(huán)境。

圖3 采用的帶有金屬偽影的體模與真實臨床CT圖像。(a)帶有金屬偽影的體模CT圖像;(b)帶有金屬偽影的臨床CT圖像Fig.3 The original images.(a)phantom image;(b)clinical CT image

非均值修補算法能有效利用圖像冗余信息，由于要遍歷整幅圖像，因此算法效率不高。其濾波效果對遍歷路徑、窗口半徑大小、衰減常數(shù)和判別閾值δ(見1.3中步驟4)十分敏感。筆者進行大量實驗，最終采用Rsim=5，h=2.5對仿真體模圖像進行修補，采用Rsim=3，h=1.2對臨床圖像進行修補。但實驗中由于相似程度過低，而使用線性插值的情況也有出現(xiàn)(見1.3中步驟4)。判斷標(biāo)準(zhǔn)為：每次固定目標(biāo)窗，遍歷所有相似窗而求得的最大權(quán)重wmax，當(dāng)wmax≥δ=0.3時采用非局部均值修補，否則采用雙線性插值修補。

2 結(jié)果

分別采用本算法和文獻[2]中算法校正結(jié)果如圖4所示。其中的ASF濾波強度主要由此函數(shù)中的濾波掩模大小、NAST中的高斯卷積核平滑強度(濾波尺度)決定。使用較大的掩模，濾波強度大，計算較慢。高斯核的平滑尺度越大，平滑效果越好。圖4中(a)和(c)采用文獻[2]中算法，作為對比。圖4中(b)和(d)選用掩模半徑為3，平滑尺度為1.5的ASF濾波。高強度ASF濾波可以有效去除條狀偽影，但過大的強度使得復(fù)雜的臨床圖像中被偽影污染的肌肉、軟組織等一律被濾除，因此必須確定合適的濾波強度。所有結(jié)果圖像與原始待處理的圖像采用相同的窗寬。

從圖4中(a)和(c)對于臨床圖像的校正可以看出，文獻[2]中算法可以在一定程度上濾除條狀偽影，但金屬邊緣處仍有殘余。對比之下，本校正算法可以在有效消除條狀偽影前提下，最大限度使原有的圖像不被破壞，從而保證校正后的圖像仍具有臨床價值。對比圖4(a)、(b)，對于水的重建，本算法結(jié)果(圖4(b))勻質(zhì)區(qū)域的一致性更高。對比圖4(c)、(d)，對于重建圖像中的條狀偽影的抑制更好。

圖4 分別采用文獻[2]和本算法對體模、臨床圖像的偽影消除后的結(jié)果。文獻[2]算法用于體模圖像;(b)本算法對用于體模圖像;(c)文獻[2]算法用于臨床圖像;(d)本算法對用于臨床圖像Fig.4 The corrected phantom and clinical CT image using the algorithm in Ref[2]and our proposed algorithm.(a) the corrected phantom using the algorithm in Ref[2];(b)the corrected phantom using our proposed algorithm;(c) the corrected clinical image using the algorithm in Ref[2];(d) the corrected clinical image using our proposed algorithm

3 討論和結(jié)論

所提出的算法也存在不足。如何有效區(qū)分圖像邊緣和條狀偽影是本算法亟待解決的問題。另外，非局部均值修補算法雖然可以有效利用圖像的全局信息，但算法效率不高(相比插值等方法)。進一步的工作希望構(gòu)造修補誤差泛函，通過求解其極值尋找最佳的窗口大小和衰減常數(shù)，使得非局部均值修補達到更好的效果。

CT圖像金屬偽影的消除是一個重要課題。本研究采用自適應(yīng)方向性前置濾波器對含有偽影的圖像進行處理，在一定程度上抑制條狀偽影。其后采用均值漂移分割出圖像中的金屬成分。最大互信息熵分割算法能夠較為有效地提取金屬物周圍偽影。在獲取原始圖像和金屬成分的弦空間數(shù)據(jù)后，本算法采用非局部均值結(jié)合線性插值方法對弦空間圖像進行修補。當(dāng)相似性低于給定閾值時，直接采用線性插值方法由待修補區(qū)域周圍數(shù)據(jù)對其進行修補。實驗表明，這樣的結(jié)合有助于減少弦空間修補時的誤差，最后采用濾波反投影得到校正后的圖像。

[1]Manzke R， Grass M， Hawkes D. Artifact analysis and reconstruction improvement in helical cardiac cone beam CT[J].IEEE Trans on Medical Imaging，2004，23(9)： 1150-1164.

[2]Matthieu B and Lothar S.Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[J].Medical Physics，2006，33(8)：2852-2859.

[3]Ullrich K.The edge and junction detection with an improved structure tensor[J].Pattern Recognition，2003，36(2)：25-32.

[4]Lu Qingwen，Chen Wufan.Image segmentation based on mutual information[J].Chinese Journal of Computers.2006，29(2)：296-301.

[5]Wong A，Orchard J.A nonlocal means approach to exemplarbased inpainting[J].ICIP，2008，14(3)： 2600-2603.

[6]Buades A，Coll B，Morel JM.A review of image denoising algorithms，with a new one[J].Multiscale Modeling and Simulation，2005，4(2)：490-530.

[7]Buades A.Image and film denoising by non-local means[D].Illes Balears： University of Les Illes Balears，2006.

[8]Wang Wei， Gao Jinghuai， Li Kang. Structure-adaptive anisotropic filter with local structure tensors[C]//IEEE Computer Society. Second International Symposium on Intelligent Technology Application： Shanghai： IEEE，2008：1005-1010.

[9]Brox T，Weickert J，Burgeth B，et al.Nonlinear structure tensors[J].Image and Vision Computing.2006，24：41-55.