基于頂帽運算和各向異性擴散的CT圖像金屬偽影去除

牛曉偉, 孔慧華, 邸云霞

(中北大學 數(shù)學學院 山西省信息探測與處理重點實驗室, 山西 太原 030051)

0 引言

自1895年倫琴發(fā)現(xiàn)X射線以來,計算機斷層掃描技術(shù)(computed tomography,CT)已廣泛用于臨床診斷.現(xiàn)代CT掃描儀可以提供具有可靠解剖信息的高分辨率圖像,然而,如果掃描場中存在金屬植入物(如牙科填充物或修復體),金屬偽影往往會影響CT圖像質(zhì)量并降低其診斷價值[1].導致金屬偽影產(chǎn)生的原因可歸納為射束硬化、部分容積效應(yīng)、噪聲、散射和運動等[2],這些問題會導致檢測器接收到的投影數(shù)據(jù)損壞或丟失.在圖像重建中,損壞的數(shù)據(jù)表現(xiàn)為或亮或暗的偽影,除此之外,金屬植入物的大小、形狀、位置也會影響金屬偽影的情況,這會降低CT圖像的質(zhì)量,影響金屬植入物周圍組織結(jié)構(gòu)的診斷信息.

近年來,越來越多的金屬偽影去除(metal artifact reduction,MAR)方法用來校正金屬偽影,它們大致可分為迭代重建、正弦圖校正和混合方法[3].迭代重建和混合方法在抑制金屬偽影方面表現(xiàn)出良好的性能,其中最大似然和代數(shù)重建技術(shù)是兩種最常見的方法.如Kubo等[4]比較了基于模型的迭代重建和混合迭代重建與金屬偽影去除算法相結(jié)合時改善牙齒金屬偽影CT圖像的能力,發(fā)現(xiàn)混合迭代重建算法顯著減少了頸部CT圖像的金屬偽影.但由于迭代算法在臨床應(yīng)用中受到限制,計算成本高,耗時嚴重等問題,因此,在臨床實踐中應(yīng)用有限.

在正弦圖校正方法中,大致可分為兩組:基于插值的方法和基于先驗的方法.基于插值的方法使用插值技術(shù)來估計金屬投影區(qū)域的替代數(shù)據(jù),例如多項式插值、小波插值或線性插值等[5-7].線性插值(linear interpolation,LI)是最常用的方法,Kalender等[6]在1987年首次提出了關(guān)于線性插值的金屬偽影去除算法,使用金屬軌跡兩側(cè)未損壞的數(shù)據(jù)對金屬投影區(qū)域進行線性插值,這樣往往會造成組織信息的丟失,而且可能會產(chǎn)生嚴重的二次偽影.為了克服基于插值方法的缺點,基于先驗的方法使用來自先驗圖像的前向投影來估計替代數(shù)據(jù),先驗圖像的結(jié)構(gòu)信息可以添加到金屬投影區(qū)域,這樣可以獲得更好的估計[8-15].Wang等[8]提出了一種基于融合先驗的MAR(fusion prior-based MAR,FP-MAR),將濾波后的線性插值圖像和去除金屬部分的原始重建圖像融合生成一個先驗圖像,保留了更多的組織信息,但會殘留白色偽影.Meyer等[9]基于初始圖像的多閾值分割獲得的三值圖像的正向投影來執(zhí)行歸一化,提出了歸一化金屬偽影去除(normalized metal artifact reduction,NMAR)算法,以便更好地保持空氣和水等物體之間的對比度,對于帶有金屬植入物的骨骼結(jié)構(gòu)有相當大的改進.Yakdiel等[12]利用Mumford-Shah圖像模型和L0梯度最小化方法,通過高斯濾波、Parisotto和Schoenlieb修復方法三種變換獲得先驗圖像,可以有效減少條紋金屬偽影,避免金屬植入物周圍出現(xiàn)新偽影.

由于一般基于先驗的方法通常會在最終校正圖像中造成偽影殘留和金屬附近骨骼結(jié)構(gòu)丟失,如基于線性插值的NMAR算法等,無法取得理想的校正效果.因此,本文提出一種基于頂帽運算和各向異性擴散[16]的歸一化金屬偽影去除算法(top hat and anisotropic diffusion-NMAR,THAD-NMAR).該算法對CT圖像執(zhí)行頂帽運算并與線性插值圖像進行融合,可以保留大部分骨骼信息和細節(jié)特征,且各向異性擴散濾波[16]具有保留圖像邊緣細節(jié)同時減少噪聲的作用,由此可以獲得質(zhì)量更好的先驗圖像.臨床圖像測試結(jié)果表明,與現(xiàn)有的常規(guī)圖像域金屬偽影去除算法相比,本文所提出的THAD-NMAR方法在抑制偽影和組織特征保留方面表現(xiàn)更好.

1 理論與方法

含金屬物體的重建圖像常常表現(xiàn)為強烈的黑色或白色金屬偽影,或像骨骼一樣明亮,或像空氣一樣暗,經(jīng)常會出現(xiàn)分割錯誤的情況,嚴重的偽影會覆蓋感興趣區(qū)域,干擾正確診斷.而先驗圖像可以在一定程度上彌補常規(guī)算法的不足,因此先驗圖像的質(zhì)量在基于先驗的MAR中至關(guān)重要,為了克服由線性插值引起的偽影殘留和組織結(jié)構(gòu)丟失的問題,本文使用了一種數(shù)學形態(tài)學算法.主要思想是通過對CT圖像執(zhí)行頂帽運算來恢復因線性插值方法丟失掉的骨骼信息,然后使用各向異性擴散濾波去除殘留偽影,從而得到質(zhì)量更好的先驗圖像.本文算法主要分為以下四個步驟.

1.1 金屬分割

當被掃描物體中包含多個金屬材料時,會導致重建的CT圖像中出現(xiàn)嚴重的金屬偽影,如圖1(a)所示,黃色框為黑色偽影,紅色框為白色偽影.首先需要預先從圖像中分割出這些金屬來得到金屬軌跡,而閾值分割法具有速度快的優(yōu)點,能準確識別出金屬部分,可以滿足大多數(shù)CT圖像的金屬分割.由于CT圖像中不同的組織具有不同的CT值,例如空氣為-1 000 HU左右、脂肪為-90～-70 HU、水為0 HU左右、軟組織為20～40 HU,骨骼為300～1 000 HU,而金屬的CT值遠遠大于2 000 HU,因此將分割的閾值選擇為2 500 HU,可以得到CT圖像的金屬部分,對金屬部分前向投影得到金屬軌跡.

圖1 黑色偽影校正圖像

1.2 分類線性插值

在投影域中對原始圖像投影所有角度下的金屬軌跡部分使用線性插值,公式為:

(1)

式(1)中:gβ為角度β下校正后的投影值,[mβ,nβ]為該角度下的金屬投影所在的區(qū)間,即插值區(qū)間,坐標α為對應(yīng)的插值點,且α∈[mβ,nβ].然后對線性插值的投影重建得到線性插值圖像Ili,線性插值圖像雖然可以有效去除白色和黑色金屬偽影,但會引入嚴重的二次偽影,造成金屬植入物骨骼丟失的情況,如圖1(b)所示.通過下面公式融合可以得到無黑色偽影圖像:

Ino_black=Iorign+Iblack

(2)

其中,

(3)

式(3)中:Iblack是黑色偽影的補償信號,Iorign是原始圖像,通過式(2)可以得到無黑色偽影圖像Ino-black.無黑色偽影圖像通過添加補償信號只是去除了原始圖像的黑色偽影,并沒有校正白色偽影,可以保留原始圖像中的所有細節(jié)特征,如圖1(d)所示.

1.3 先驗圖像生成

1.3.1 基于頂帽運算的骨骼修復

MAR中的插值算法常常會導致最終校正圖像中金屬附近的組織結(jié)構(gòu)丟失.數(shù)學形態(tài)學是圖像處理過程中一種非常重要的方法,可以提取圖像中的分量信息,包括膨脹、腐蝕、開運算和閉運算四種基本運算,開運算可以去除相對于結(jié)構(gòu)元素較小的明亮細節(jié),保持整體的灰度和較大的明亮區(qū)域不變,頂帽運算正是用來提取這些被移除的部分,選擇不同的結(jié)構(gòu)元素會導致不同的分割結(jié)果,即提取出不同的特征.對無黑色偽影圖像Ino-black執(zhí)行頂帽運算,可以恢復丟失掉的骨結(jié)構(gòu).頂帽運算的定義為:

假設(shè)原始灰度圖像為I,使用半徑大小為d的圓盤結(jié)構(gòu)元素z進行頂帽運算,公式如下:

Itophat=I-(I°z)

(4)

其中,

I°z=(I?z)⊕z

(5)

式(4)中:Itophat是頂帽運算提取到的信息,?表示腐蝕,⊕表示膨脹,結(jié)構(gòu)元素z的大小d依據(jù)骨骼信息的情況而定.式(5)為開運算,用°表示,開運算是對圖像先腐蝕,再膨脹.腐蝕的公式為:

I?z=min[I(x+i,y+j)-z(i,j)]

(6)

式(6)中:(x+i,y+j)是定義在目標I上的坐標,(i,j)是定義在結(jié)構(gòu)元素z上的坐標.膨脹的公式為:

I⊕z=max[I(x+i,y+j)+z(i,j)]

(7)

先驗圖像的質(zhì)量在基于先驗的MAR中起著至關(guān)重要的作用,而傳統(tǒng)的插值方法往往會由于分割錯誤,在插值過程中丟失骨骼等組織結(jié)構(gòu),從而影響到圖像質(zhì)量.通過頂帽運算可以恢復丟失的骨骼信息,從而獲得更好的先驗圖像.

將無黑色偽影圖像減去開運算的結(jié)果,可以提取到圖像的骨骼部分,這是由于偽影的CT值一般小于300 HU,而骨骼的CT值為300～1 000 HU,因此選擇閾值200 HU可以去除其他細節(jié)(包含噪聲).然后對線性插值圖像和無黑色偽影圖像進行融合,融合過程為:

(8)

式(8)中:Itophat表示CT圖像頂帽運算的結(jié)果,如圖2(a)所示,Ifusion表示融合后的圖像,它會將線性插值丟失的骨骼信息恢復,同時去除圖像中的大部分偽影,如圖2(b)所示.

圖2 頂帽運算骨骼恢復圖像

1.3.2 各向異性擴散濾波

融合后的CT圖像往往會殘留一部分偽影,因此需要進一步處理.傳統(tǒng)的線性濾波算法在去噪時通常會犧牲掉圖像的邊緣信息,而基于非線性偏微分方程的各向異性擴散方法會彌補這個缺陷,具有保留圖像邊緣細節(jié)同時減少噪聲的作用,可以增強圖像邊緣的對比度.Perona和Malik給出的各向異性擴散方程(P-M方程)的離散表達式為[16]:

(9)

(10)

對融合后的圖像Ifusion執(zhí)行各向異性擴散濾波,可以得到先驗圖像Iprior,如圖3所示.

圖3 各向異性擴散濾波圖像

1.4 歸一化金屬偽影去除

基于原始CT圖像、先驗圖像和金屬部分的投影,執(zhí)行Meyer等[9]提出的歸一化和反歸一化操作,對受金屬影響正弦圖進行校正,過程如下:

(11)

式(11)中:M(·)表示投影的線性插值,R表示拉東變換,Pcorrect表示最終校正的正弦圖.最后通過濾波反投影算法重建得到最終的校正圖像Icorrect.

2 實驗結(jié)果與分析

2.1 實驗設(shè)置

為了驗證本文算法在去除金屬偽影方面的有效性,使用臨床CT圖像進行測試.本文選擇了6個臨床病例來測試所提出的THAD-NMAR算法在實際臨床中應(yīng)對不同種類金屬偽影的能力.所有原始未校正金屬偽影圖像均來自REVISION RADIOLOGY(http://www.revisionrads.com/cgi-bin/MDT_cgi).圖4中包含了三組金屬結(jié)構(gòu)簡單的圖像,圖5包含了三組金屬結(jié)構(gòu)較為復雜、偽影嚴重的圖像.本文選取線性插值法、NMAR、LI-NMAR、金屬刪除技術(shù)(metal deletion techniq-ue,MDT)和FP-MAR算法作為比較算法,其中FP-MAR算法的參數(shù)均設(shè)置為原文中的值,本文算法各向異性擴散濾波的參數(shù)都是根據(jù)經(jīng)驗選擇,可以達到最優(yōu)效果,其中迭代次數(shù)t設(shè)置為100,熱傳導系數(shù)K設(shè)置為15,控制擴散總體強度的λ設(shè)置為0.25,可滿足所有的實驗需求.圖像大小均為512×512,CT圖像窗位和窗寬均設(shè)為C=500/W=3 000 HU.

2.2 實驗分析

圖4顯示的是臨床數(shù)據(jù)金屬植入物對比方法的測試結(jié)果,案例1是腦動脈瘤植入物,案例2和案例3是深部腦刺激器,且參數(shù)d的值均為0.圖4是幾個金屬偽影較輕的情況,所有的方法都可以消除大部分偽影,其中線性插值法在校正后會引入二次偽影,NMAR算法會因為偽影較嚴重而分類錯誤最終殘留下一部分黑色偽影,LI-NMAR算法較NMAR算法校正效果會好很多.而FP-MAR算法由于是通過融合濾波后的線性插值圖像和去金屬的原始偽影圖像來恢復丟失掉的骨信息,因此會殘留下白色偽影,如圖4(f)中的案例3所示.本文的THAD-NMAR算法通過有效的偽影抑制和骨骼結(jié)構(gòu)保留獲得了最佳的圖像質(zhì)量.

圖4 臨床數(shù)據(jù)金屬偽影去除方法結(jié)果對比

圖5是更為復雜的情況,包含了多個不規(guī)則或大型金屬,其中案例4是胸骨固定螺釘,案例5是單側(cè)髖關(guān)節(jié)假體,案例6是椎弓根螺釘,d值分別為10、12、35.線性插值法、NMAR算法和LI-NMAR算法雖然去除了大部分的金屬偽影,也損失了金屬周圍的骨骼信息,如圖5(b)、(c)、(d)所示,這樣的結(jié)果不利于診斷.MDT方法有很好的偽影抑制效果,但也會缺少一部分金屬附近的組織結(jié)果并殘留白色偽影,如圖5(e)中的案例4.FP-MAR算法在很大程度上抑制了金屬偽影和保留大部分的骨骼信息,卻會殘留下一部分明顯的白色偽影,如圖5(f)中的案例4和案例6所示,因此具有一定的局限性.顯然,本文提出的THAD-NMAR算法在抑制偽影方面優(yōu)于其它MAR算法,在恢復那些靠近金屬物體的組織結(jié)構(gòu)方面也比MDT和FP-MAR算法表現(xiàn)更好.

圖5 臨床數(shù)據(jù)金屬偽影去除方法結(jié)果對比

金屬偽影或像骨骼一樣明亮,或像空氣一樣暗,就經(jīng)常會出現(xiàn)分割錯誤的情況,最終會殘留下一部分黑白偽影,因此擁有更好的先驗圖像顯得尤為重要.LI-NMAR將線性插值圖像作為先驗圖像,由于在金屬較多時線性插值圖像會造成明顯的殘留偽影和骨骼缺失,最終的校正圖像也會出現(xiàn)骨骼缺失和偽影殘留的情況.FP-MAR算法也有一定的局限性,例如它可能會將金屬物體周圍的一些明亮偽影錯誤地合并到先驗圖像中作為骨骼結(jié)構(gòu),因此可能導致最終校正后殘留下白色偽影.而本文的THAD-NMAR算法恰好彌補了這一缺陷,它可以準確地提取出原始圖像中的骨骼信息.

由于圖4中的MAR算法不會出現(xiàn)骨骼缺失的情況,因此參數(shù)d的值均設(shè)為零,而圖5中的CT圖像有較為嚴重的金屬偽影,所以在常規(guī)MAR算法校正后會出現(xiàn)偽影殘留和骨骼信息丟失的情況,依據(jù)不同的情況需要設(shè)置不同的d值.圖6是從無黑色偽影圖像中提取到的骨骼信息,可以看到,不同的d值會提取到不同的信息,當d值為10和15時可以提取到程度相同的骨骼和金屬部分.d值過小,提取到的骨骼信息可能會不完整,造成組織信息丟失;d值過大,會提取到原始CT圖像中的白色偽影和其他信息,進而造成偽影殘留.因此需要手動設(shè)置選取合適的d值,當骨骼少或較細時d值應(yīng)小一點,當骨骼較多或復雜時d值要大一點.

圖6 案例4頂帽運算不同d值提取到的骨骼信息

圖7和圖8顯示了沿案例4原始CT圖像中第212行(紅色線條)和296行(黃色線條)得到的各種算法校正圖像的強度分布曲線,用來對比黑色偽影和白色偽影的校正情況.從圖7可以看到,當前黑色偽影較為嚴重,每種算法都可以較好的去除黑色偽影,保持正常的CT值.

從圖8可以看到,當前白色偽影較為嚴重,CT值遠遠大于正常值,除FP-MAR算法外都可以很好地校正骨骼附近的白色偽影,但LI算法、NMAR算法、LI-NMAR算法會造成骨骼丟失的情況,如圖8(a)中的紅色虛線框所示,而THAD-NMAR算法不僅將CT值恢復到了正常的CT值附近,而且可以保留骨骼等重要信息,因此本文算法在去除白色偽影和保留骨骼信息方面要優(yōu)于其他算法.

圖7 案例4第286行校正CT圖像的強度分布

圖8 案例4第212行校正CT圖像的強度分布

3 結(jié)論

CT圖像中的金屬偽影一直是影響CT圖像質(zhì)量的關(guān)鍵因素之一,這些金屬偽影不僅模糊了CT圖像,而且影響了疾病的診斷效果,甚至影響了CT成像系統(tǒng)的進一步發(fā)展.而一般基于先驗的金屬偽影去除方法通常會在最終校正圖像中造成偽影殘留和金屬附近骨骼結(jié)構(gòu)丟失,為了解決這個問題,本文提出了一種基于頂帽運算和各向異性擴散濾波的歸一化金屬偽影去除算法.

通過頂帽運算來恢復因線性插值引起的組織信息丟失問題,然后通過各向異性擴散濾波可以得到先驗圖像,最后基于原始CT圖像投影、先驗圖像投影和金屬軌跡執(zhí)行歸一化和反歸一化操作得到最終的校正圖像.頂帽運算通過調(diào)整參數(shù)d可以恢復丟失掉的骨結(jié)構(gòu),而各向異性擴散濾波具有保留圖像邊緣細節(jié)同時減少噪聲的作用.

與其他經(jīng)典的MAR算法相比,在存在多個或較為復雜的金屬植入物時,本文算法可以去除CT圖像中強烈的暗和亮的金屬偽影,不會引入新的偽影,并且可以完整地保留金屬周圍的組織結(jié)構(gòu),適合較大金屬植入物周圍存在骨骼等結(jié)構(gòu)的情況.