骨邊界增強濾波的圖割算法

石志良,范偉楠,甘梓博,袁 瓊

(1.武漢理工大學機電工程學院,湖北 武漢 430070)(2.湖北第二師范學院計算機學院,湖北 武漢 430205)

醫學圖像具有灰度廣、對比度低等特點,快速精準地分割骨組織一直是難題. 臨床通常采用閾值分割和區域生長相結合的方法進行手工分割,不僅需要大量人工交互,而且需要操作人員具備深厚的專業知識,分割結果易受主觀因素影響,具有不可重復性. 因此,研究一種快速、穩定的自動分割算法具有重要意義.

常見的骨自動分割算法分為機器學習和非機器學習兩類. 機器學習方法,可實現令人滿意的精度和自動化,但技術人員需要對特定器官進行大量標注和反復訓練,過程復雜、耗時,且僅適用于單個器官. 非機器學習方法指借助圖像中固有的灰度、梯度、紋理等信息,根據一定的相似性準則,將圖像分割成各具特性的區域. 主要方法包括閾值分割[1]、區域生長[2-3]、水平集[4-5]、聚類[6-8]、圖割[9-20]等,這類算法具有通用性的優點,適用于更廣泛的分割問題,因此研究一種非機器學習的自動分割算法具有不可替代的意義.

近年來,圖割作為一種非機器學習的方法,具有穩定、全局最優、可擴展性好等優點,在圖像分割領域占據著重要的地位. 其中最為經典的是Boykov等[10-11]提出的基于圖像強度的交互式圖割算法. 然而,由于算法僅考慮圖像強度,在分割骨組織時魯棒性較差,并且采用交互式的方法增加圖割約束,效率低,無法滿足自動分割的需求. 為提高算法的魯棒性和效率,許多研究人員對圖割算法進行改進,先后提出了基于K聚類的圖割[9]、基于雙通道的圖割[12]、基于形狀約束的圖割[13]、基于分水嶺技術的交互式圖割[14]、基于顯著性的圖割[15]、基于特征檢測的圖割[16]、基于圖像數據映射的圖割[17]、結合U-net網絡的圖割[18]、金字塔圖割[19]、基于區域生長的圖割[20]等相關算法.

雖然上述算法在一定程度上提高了圖割的魯棒性,但未能徹底解決算法的交互問題,尤其在分割對比度低的骨邊界時效果較差. 因此,研究提出一個基于Hessian矩陣的骨骼邊界增強濾波器,在此基礎上建立新的能量函數,求解能量函數得到初步分割結果,并對初步結果進行后處理,實現骨組織的自動分割.

1 傳統圖割算法

圖割算法的基本思想是將N維圖像的分割問題轉換為求解最小能量函數問題. 能量函數的解與圖模型中的割線唯一對應,并將圖分為目標和背景兩部分. Boykov等[10]提出的能量函數為

E(A)=λ·R(A)+B(A).

(1)

式中,系數表示R(A)相對于B(A)的重要程度,A=(A1,…,Ap,…A|P|)是一組二進制向量,P是圖中所有像素的集合,Ap對應P中像素p,為目標點(表示為1)或背景點(表示為0),當所有Ap已知時,向量A即為圖的一種分割方式.

R(A)為區域項,定義為

(2)

式中,RP(AP)表示將像素標記為目標或背景的懲罰,分別用RP(“obj”)和RP(“bkg”)表示.

B(A)為邊界項,定義為

(3)

式中,當Ap=Aq時,δ(Ap,Aq)=1,否則δ(Ap,Aq)=0;B{p,q}表示相鄰像素p,q的一致性懲罰,這種一致性可以基于圖像強度、梯度或其他標準.當像素p,q相似(如在圖像強度上相近)時B{p,q}很大,反之接近于0.N表示由像素p及其相鄰像素點構成的鄰域系統.

為獲得全局最優的結果,用戶需要標記目標種子點和背景種子點以增加圖模型的輸入.通常,將區域項和邊界項稱之為軟約束,種子點稱之為硬約束.

圖1為傳統圖割算法示意圖,圖1(a)為灰度圖像,建立的加權無向圖G=(V,E),如圖1(b)所示,V為圖中節點,表示為

圖1 圖割算法示意圖

V=P∪{S,T},

(4)

式中,P為像素點集合;S,T為終端節點,分別稱之為源、匯.

E為圖中邊,表示為

E={p,q}∪ {p,S}∪{p,T},

(5)

式中,p,q為像素點;{p,q} 表示相鄰像素構成的邊,稱為n-links;{p,S} 表示像素點與源構成的邊;{p,T} 表示像素點與匯構成的邊.{p,S} ,{p,T} 統稱為t-links.

邊的權由軟約束和硬約束共同決定,{p,q} 的權為B{p,q};{p,S} 的權分為3種情況,當p為背景種子點時其值為0,p為目標種子點時其值為1+max(B{p,q}),p為非種子像素點時其值為λ·Rp(“bkg”);類似的{p,S} 的權為,當p為背景種子點時其值為1+max(B{p,q}),p為目標種子點時其值為0,p為非種子像素點時其值為λ·Rp(“obj”).

根據上述原理,用最大流算法求解能量函數,使得t-links和n-links權值最小. 如圖1(c)所示,曲線1即為所求割線,最終分割結果如圖1(d)所示.

2 本文方法

2.1 骨邊界增強濾波的圖割模型

算法總體流程如圖2所示,主要由骨邊界增強濾波、能量函數建立和后處理3個部分組成.

圖2 本文算法流程圖

骨邊界增強濾波輸出結果的好壞直接決定骨邊界的清晰度,進而影響圖像分割質量,因此,濾波器的設計是本文第一個關鍵步驟. 如何將濾波結果融合到圖割約束項中,建立包含圖像強度與結構特征的能量函數,是本文第二個關鍵步驟. 此外,由于CT圖像固有的部分體積效應,分割后的骨組織可能存在粘連的情況,如何設計一種后處理方法將相鄰骨分離,提高分割的魯棒性和效率,是本文第三個關鍵步驟. 針對這3個關鍵步驟對圖割算法進行改進,以期提高圖像的分割精度與自動化程度.

2.2 骨邊界增強濾波器

骨分割難點之一是關節處骨邊界對比度低,因此僅基于圖像強度特征的分割方法,不能取得理想的分割效果. 為提高骨邊界對比度,使用經典的邊緣增強濾波器(非銳化遮蔽濾波)銳化圖像,以突出CT圖像中灰度突變的區域. 圖3表示股骨切片圖像非銳化遮蔽處理后的效果,由圖可知,骨邊界強度較原圖得到提升,但相鄰骨之間仍存在連接,對分割結果產生干擾.

圖3 非銳化遮蔽效果圖

為此,基于Hessian矩陣設計一種增強骨邊界的濾波器,稱為骨邊界增強濾波.

Hessian矩陣是由多元函數的二階偏導數構成方陣,主要用來描述函數局部曲率. 以二維圖像為例,令表示圖像像素關于坐標的函數,將其在泰勒展開得

(6)

則其Hessian矩陣為

(7)

在二維圖像中,Hessian矩陣是二維正定矩陣,具有2個特征值和2個特征向量,其中特征值表示圖像在2個特征向量所指方向上的各向異性.

利用特征向量和特征值(λ1,λ2)繪制橢圓,如圖4(a)所示.當λ1=λ2時,各向同性最強,對應圖4(b)中斑點結構;當λ1=1,λ2=0或λ1=0,λ2=1,各向異性最強,對應圖4(c)中線性結構. Hessian矩陣特征值與圖像特征對應關系見表1.

表1 二維Hessian矩陣特征值與圖像特征對應關系

圖4 二維圖像結構特征

二維圖像求二階導的一般方法為

fxx(x,y)=f(x,y)-f(x-Δx,y)-(f(x-Δx,y)-f(x+2Δx,y)),

(8)

顯然此方法只包含自身在內的3個像素信息,易受圖像局部信號干擾.根據線性尺度空間理論,對一個函數求導等于該函數與高斯函數導數的卷積,其表達如下

(9)

式中,σ表示高斯函數的標準差,稱為尺度,fxx(x,y;σ)表示在圖像尺度σ下的二階偏導數.高斯模板包含周圍矩形范圍內所有像素信息,因此能顯著降低求導誤差.

由于在一幅圖像中,同種結構具有不同大小,如粗細不同的線,故同一尺度無法適用于整張圖像.為此,采用多尺度方法,即對同一點使用不同尺度的高斯模板進行卷積,選擇各向異性最強的結果作為該點的輸出.

將上述理論推導至三維圖像,得尺度下的三維Hessian矩陣為

(10)

式中,表示三維圖像在體素v的二階偏導數.設矩陣3個特征值為λ1,λ2,λ3(|λ3|≥|λ2|≥|λ1|).以3個特征值為半軸可得特征橢球如圖5所示.

圖5 特征橢球

顯然,當圖像各向同性強時橢球將被壓縮為球體;當有兩個特征值的絕對值接近1,另一個接近0時,橢球將被壓縮為柱體;當有兩個特征值的絕對值接近0,另一個接近1時橢球將被壓縮為片狀. 因此可將三維圖像的特征結構分為片狀、管狀和斑點狀3種,如圖6所示,與特征值的對應關系見表2.

表2 三維Hessian矩陣特征值與圖像特征對應關系

圖6 理想結構

其中噪聲狀是斑點狀的特例,其特點是特征值均為0;骨邊界是由骨外層的皮質骨組成,由表可得皮質骨對應的結構為片狀,因此使用Hessian矩陣提取骨邊界并對其進行增強是可行的.

為量化CT圖像體素v局部區域與片狀結構的相似性,可定義圖像在尺度下的評估函數

(11)

式中,α,β,γ為常數,R1,R2,R3表示與片狀、管狀、斑點狀、噪聲狀有關的項,稱之為測量因子,分別定義為

(12)

根據定義計算出不同結構對應的測量因子R1,R2,R3的值,見表3.

表3 圖像結構與測量因子的對應關系

如上文所述,同一尺度無法適用于整張圖像,為提高評分函數的魯棒性,計算不同尺度下的評分最大絕對值,并將最大特征值納入函數中,得

(13)

式中,∑為不同尺度構成的集合.

圖7為股骨CT切片經骨邊界增強濾波器處理效果圖,由圖7(b)可知,處理后相連的股骨與髖骨分離開來,骨邊界清晰度度提高,證明本文提出的基于Hessian矩陣的骨邊界增強濾波器是有效的.

圖7 濾波處理效果圖

2.3 約束項計算

2.3.1 改進區域項函數

傳統的區域項函數僅基于圖像強度,由于骨邊界薄弱、間隙狹窄、骨小梁處強度低,不能取得很好的分割結果.為此考慮結合圖像強度與評分函數,定義骨區域和背景區域.

在CT圖像中圖像的強度經過線性變換可以得到CT值,其單位為Hounsfield Unit(HU),CT值與圖像組織的對應關系,見表4.

表4 不同組織對應HU值

由表4可得,骨組織的CT值一般在400HU以上,軟組織的CT值一般不超過150HU,空氣和脂肪的CT值一般小于-80HU,結合評分函數,定義估算的骨組織區域Ebone和背景區域Ebkg分別為

Ebone={x∈Ω|I(x)≥400且S(x)>0},Ebkg=lcc({x∈Ω|I(x)<-80}),

(14)

式中,表示體素的CT值,表示高強度的皮質骨,表示為排除骨間通道中的軟組織體素,lcc表示輸入二進制參數的最大連通分量,{x∈Ω|I(x)<-80}表示脂肪和空氣. 根據式(14)對股骨切片圖像進行處理,結果如圖8所示,其中圖8(b)為非骨部分的像素區域記為Ebkg,圖8(c)為非骨部分的像素區域記為Ebone.

圖8 區域估計

由區域項定義可知,區域項函數中表示將體素p標記為目標點(或背景點)的懲罰.分割骨組織時,在體素p是骨骼或在體素p是背景且時懲罰小,其他情況懲罰高.結合估算的區域,建立新的區域項函數為

(15)

新的區域項函數綜合評分函數值和圖像強度,相較于單一的基于圖像強度的區域項函數具有更明顯的特征和更高的可靠性.

2.3.2 改進邊界項函數

同傳統區域項函數一樣,僅基于圖像強度的邊界項函數不能處理對比度低的骨邊界,因此考慮使用評分函數設計新的邊界項函數.

邊界項中表示相鄰體素p,q的一致性懲罰,p,q越相似(如在圖像強度上相近)懲罰越大,反之則小.因此,基于評分函數提出新的邊界一致性標準,當評分函數值相近時懲罰大,反之則小.定義邊界項函數為

(16)

根據式(1)、(2)、(3)、(15)、(16),得能量函數E(A)為

(17)

根據上述原理對CT圖建立圖模型,使用最大流算法求解能量函數得出圖模型中對應割線,該割線將圖像分割為前景背景兩部分,前景即為所求骨組織. 圖9為分割股骨切片的示意圖,圖9(a)為原始切片;圖9(b)為圖割后的二值結果,圖中可見,切片被分為前景、背景兩部分,其中前景由兩個不連通區域A和B構成,分別代表髖骨和股骨,兩者邊界清晰、互不交叉;圖9(c)為將原圖與二值圖相疊加的效果圖,圖中可見,骨組織與軟組成邊界清晰,骨分割質量好,證明本文提出的基于骨邊界增強濾波的圖割算法是可行的.

圖9 股骨切片初步分割結果

2.4 后處理

雖然上述骨邊界增強濾波能夠提高骨邊界清晰度,但由于CT圖像固有的部分體積效應,仍會導致初步分割結果中出現相鄰骨組織粘連的情況,如圖10所示,圖10(a)是使用上文方法獲得的初步二值結果,圖中可見股骨區域與髖骨存在細微的連接點,導致兩骨無法徹底分離.為分離相鄰的骨骼,提出基于形態學和簡易圖割的后處理方法,將半徑為R的球形元素作為卷積核對分割后的二值圖像進行形態學腐蝕處理,如10(b)所示,對原圖進行腐蝕操作,區域C被分割為D、E和孤島3個區域.

圖10 相鄰骨分離示意圖

對于骨的分離操作,目標是在邊界體素最少時,將圖像分割為F、G兩個區域,其中,對圖像做簡易圖割處理,步驟為

(1)添加硬約束,令腐蝕操作后的區域D和F分別為圖割中的前景種子點和背景種子點;

(2)建立區域項函數,令區域項函數為

(18)

(3)為得到包含體素最少的邊界,令邊界項函數為;

(4)建立能量函數并使用最大流算法求解.

(5)將分離后的相鄰骨賦予不同CT值作為標簽.

簡易圖割處理結果如圖10(c)所示,其中圖10(b)中的孤島區域被納入區域D中,符合預期,證明本文提出的基于形態學與簡易圖割的后處理方法是有效的.

3 實驗及結果

為驗證本文算法,可對骨組織進行自動精確的分割,在CT影像數據集上進行實驗,并與閾值分割、區域生長、傳統圖割進行對比.

3.1 實驗數據集

實驗數據集分為兩組:第一組為內部數據集,來源于武漢大學中南醫院,包含30例下肢CT圖像,圖像層內間距0.78 mm,層間間距1.25 mm,圖像強度單位為HU. 以醫學專家進行人工分割結果作為金標準;第二組為公開數據集,來自TCIA,名為CT Lymph Nodes,包含176例腹部CT圖像.

3.2 實驗環境及參數設置

實驗基于Windows10系統,采用Visual Studio 2015開發平臺,配置醫學圖像處理工具包ITK、可視化工具包VTK、QT框架和最大流庫,使用C++語言開發,實現圖切割的高效計算. 實驗PC硬件配置為Intel(R)Core(TM)i7-8700K CPU@3.70GHz,16GB內存.

經過大量的實驗得出本文算法及對比算法合適參數見表5.

表5 本文算法及對比算法主要參數

3.3 評價指標

醫學圖像分割領域中評價指標眾多,能夠從不同角度對算法的預測結果進行評估.本文采用Dice系數(Dice)、精確率(Precision)、召回率(Recall),以及F1分數(F1Score)作為評價指標.

(1)Dice系數:表示兩個樣本的重疊情況,設A為人工分割金標準的圖像像素的集合,B分割算圖像像素集合,其定義為

(19)

當A與B不相交時,Dice相似性系數為0;當A與B完全相交時,Dice相似度系數為1,即系數越大,圖像分割質量越高.

(2)精確率:表示分割的骨區域中真實骨區域所占比例,其定義為

(20)

式中,TP表示骨被預測為骨,FP表示非骨被預測為骨.

(3)召回率:表示分割骨區域中的真實骨占實際骨區域的比例,其定義為

(21)

式中,FN表示骨被預測為非骨.

(4)F1分數:是精確率和召回率的調和平均數,同時考慮預測的骨區域和真實骨區域,評估結果全面,其定義為:

(22)

3.4 實驗結果分析

3.4.1 骨分割結果對比

從數據集中選取5例,使用閾值分割、區域生長、傳統圖割、骨邊界增強濾波的圖割算法對股骨進行分割,分割結果分別如圖11所示,其中(a)行是經裁剪后的原始CT圖像,(b)行是由專業醫師手工分割的實際股骨區域,作為實驗的真值,(c)行是由閾值分割算法提取的股骨區域,(d)行是由區域生長算法提取的股骨區域,(e)行是由基于強度的傳統圖割算法提取的股骨區域,(f)行是本文算法提取的股骨區域.

圖11 閾值分割、區域生長、傳統圖割和本文算法對股骨分割結果對比

對比可得,閾值分割算法分割的效果較差,是由于閾值分割的原理是將圖像中強度在一定范圍內的像素不經處理地提取出來,如四、五列所示,分割結果中不僅存在與骨組織強度相近的非骨像素,而且存在多種組織形成的不連通區域,這是導致閾值分割算法分割效果不理想的主要原因;區域生長對不連通區域的處理優于閾值分割,如第三列所示,是因為區域生長能提取圖像連通區域,但其本質是基于圖像強度,無法排除高強度的非骨像素,因此分割效果仍不理想;傳統圖割算法在骨干區域的分割效果較好,是由于骨干區域圖像對比度高,但在關節區域由于第二節介紹的原因,無法準確提取股骨區域,如第四、五列所示;本文算法提取的股骨區域與實際股骨區域更為相近,經過比較可得本文算法分割結果明顯優于其他方法.

從數據集中選取4例,分別對股骨、髖骨、脛骨、腓骨周圍區域進行裁剪.使用4種方法分割并重建,如圖12-15所示.圖12為股骨重建對比圖,其中圖12(a)為實際股骨;圖12(b)為閾值分割重建的股骨,模型存在多個不連通區域,且不能分離股骨和髖骨;圖12(c)為區域生長重建的股骨,模型不連通區域減少,但仍不能分離股骨和髖骨;圖12(d)為傳統圖割重建的股骨,模型與區域生長模型相近,不能分離股骨和髖骨;圖12(e)為本文算法重建的股骨,模型與實際股骨在形狀和大小十分貼近實際股骨,比較可得,本文算法對股骨的分割優于其他方法.圖13為髖骨重建對比圖,其中圖13(a)為實際髖骨,分割難點在于髖骨區域骨組織較多,與髖骨鄰近的有兩側股骨及中央骶骨,比較可得,本文算法能成功分割髖骨,分割效果優于其他方法.圖14、15分別為脛骨、股骨重建對比圖,比較可得本文算法對脛骨、股骨的分割效果優于其他方法.

圖12 閾值分割、區域生長、傳統圖割和本文算法股骨分割重建對比

圖13 閾值分割、區域生長、傳統圖割和本文算法髖骨分割重建對比

圖14 閾值分割、區域生長、傳統圖割和本文算法脛骨分割重建對比

圖15 閾值分割、區域生長、傳統圖割和本文算法腓骨分割重建對比

3.4.2 骨分割結果對比

為定量評價算法的可靠性,用閾值分割、區域生長、傳統圖割算法和本文算法分別分割骨組織,并計算4種骨用4種方法在評價指標Dice、Precision、Recall和F1 Score上的均值,并與Jcgs等[21]提出的3D U-Net比較,測試數據為TCIA中名為CT Lymph Nodes的數據集,評估結果見表6.

表6 4種方法骨分割評估結果

對比實驗結果,可得閾值分割算法的精確率最低,表示實際骨區域占預測骨區域的比例較低. 區域生長算法較閾值分割算法表現較好,在精確率和F1分數上有一定的提高,但各向指標均不如本文算法. 傳統圖割算法的召回率高,精確率很低,這表明傳統圖割算法對骨區域的提取不夠準確. 本文算法在保證高召回率的情況下精確率能達99.07%,對召回率和精確率的平衡較好,在綜合評價指標Dice和F1 Score上有明顯提升. 綜合4種評估結果,本文算法的整體表現優于其他方法,且與3D U-Net算法在Dice指標上相近.

4 結論

針對骨組織的自動分割進行研究,基于Hessian矩陣提出骨邊界增強濾波,增強骨邊界清晰度;將濾波結果融入圖割的能量函數中,提高圖割算法的魯棒性和效率;融合形態學和圖割提出一種后處理方法,實現相鄰骨的自動分離. 在實際應用中,該方法對骨分割處理較好,但對于其他組織的分割適應性差,因此在下一步的研究工作中,將面向多樣的人體組織進行分割研究,如血管、肺、結核等,進一步提升算法的通用性.