隨機(jī)森林和水平集的椎骨CT分割方法

陶永鵬,景 雨,頊 聰

(大連外國(guó)語(yǔ)大學(xué) 軟件學(xué)院,遼寧 大連 116044)

1 概 述

脊柱,是人體中央負(fù)重骨骼,是人體上半身結(jié)構(gòu)的中軸線(xiàn).新型掃描技術(shù)中CT是評(píng)估三維空間椎骨形態(tài)最精確的方式.椎骨CT圖像分割是診斷的基本步驟和任務(wù),如識(shí)別脊椎異常診斷[1],生物力學(xué)建?；驁D像引導(dǎo)脊柱干預(yù)[2].為了保證分割的準(zhǔn)確性,圖像引導(dǎo)脊柱干預(yù)等通常需要亞毫米級(jí)精度.但由于椎骨形狀的復(fù)雜性和結(jié)構(gòu)的變化性等特點(diǎn),傳統(tǒng)的圖像分割方法并不適用于椎骨的分割.

近年來(lái),眾多學(xué)者已提出部分針對(duì)椎骨圖像的分割算法.一些無(wú)監(jiān)督的圖像處理方法也被應(yīng)用于椎骨的分割,流域[3]和圖形[4]等處理復(fù)雜拓?fù)浣Y(jié)構(gòu)的水平集方法被應(yīng)用于椎骨分割.Lim等[5]將Willmore流程納入水平集框架來(lái)擬合表面模型的演變.黃等[6]結(jié)合邊緣和基于區(qū)域的水平集函數(shù)實(shí)現(xiàn)CT圖像中椎體的分割.李等[7]提出了一種自動(dòng)初始化水平集方法,實(shí)現(xiàn)基于混合形態(tài)濾波和高斯混合模型處理拓?fù)渥兓?Klinder等[8]提出了集成檢測(cè)框架,基于脊柱曲線(xiàn)提取和統(tǒng)計(jì)形狀模型識(shí)別和分割椎骨.Roberts等[9]使用主動(dòng)形狀模型將椎骨分成幾部分進(jìn)行分割協(xié)作,可應(yīng)用于二維射線(xiàn)圖像并可以擴(kuò)展到三維.唐利明[10]等提出了變分水平集的分割模型提高了對(duì)噪聲圖像的分割能力,但受聚類(lèi)數(shù)目的影響較大.這些模型在某種程度上達(dá)到了較高的分割準(zhǔn)確度,但初始姿態(tài)較敏感,部分并且需要對(duì)圖像手動(dòng)改動(dòng).

目前,機(jī)器學(xué)習(xí)技術(shù)已被應(yīng)用于醫(yī)學(xué)CT圖像.簡(jiǎn)[11]等使用隨機(jī)森林分類(lèi)算法成功的從存在背景干擾和光照變化的圖像中實(shí)現(xiàn)了目標(biāo)檢測(cè),但該方法并不能確定器官的范圍.Pauly等[12]應(yīng)用隨機(jī)蕨類(lèi)(random ferns),同時(shí)引入 3DLBP(Local Binary Patterns)特征,實(shí)現(xiàn)了對(duì) MR 圖像中多重解剖結(jié)構(gòu)位置和大小估計(jì).Glocker[13]采用一種有向距離圖譜對(duì)人體器官進(jìn)行檢測(cè),將有向距離圖譜作為目標(biāo)邊界形狀的隱式表達(dá),通過(guò)邊界內(nèi)外體素的正負(fù)值進(jìn)行檢測(cè).

2 本文方法

本文針對(duì)椎骨圖像分割中對(duì)初始輪廓敏感的問(wèn)題,提出基于加權(quán)隨機(jī)森林和水平集分割的椎骨CT分割方法.提取CT圖像的SIFT特征,通過(guò)隨機(jī)森林分類(lèi)回歸算法選擇距離圖譜中距離值最小的100個(gè)點(diǎn)進(jìn)行Mean-shift 聚類(lèi)分析,產(chǎn)生聚類(lèi)中心,重建操作后獲得對(duì)應(yīng)的3D距離預(yù)測(cè)圖譜,選取距離最小的點(diǎn)作為隨機(jī)森林確定的椎骨中心點(diǎn).隨后將獲得的椎骨中心點(diǎn)作為CV分割模型的初始輪廓位置,通過(guò)求解平穩(wěn)控制演化速度和噪聲敏感度的CV模型,由能量函數(shù)演化方程最小值來(lái)實(shí)現(xiàn)對(duì)椎骨CT圖像的分割.方法示意圖如圖1所示.

圖1 方法示意圖Fig.1 Schematic of the method

3 SIFT特征

SIFT特征(Scale-invariant feature transform)是一種圖像局部特征提取算法,能夠保持較好的穩(wěn)定性.在描述子的生成過(guò)程中,對(duì)于16×16的關(guān)鍵點(diǎn)鄰域,在處理梯度時(shí)進(jìn)行高斯加權(quán)處理,強(qiáng)化中心區(qū)域,淡化邊緣區(qū)域的影響.在每個(gè)4×4 子區(qū)域上計(jì)算8個(gè)方向的梯度方向直方圖,即形成一個(gè)種子點(diǎn),每個(gè)關(guān)鍵點(diǎn)描述子由4×4=16個(gè)種子點(diǎn)組成,每個(gè)種子點(diǎn)有8個(gè)梯度方向信息,形成4×4×8=128維的SIFT特征向量.

4 隨機(jī)森林

本文方法將訓(xùn)練數(shù)據(jù)點(diǎn)定義為Dk=(Sk,Ck),其中Sk是基于體素k得到的SIFT特征向量,Ck是體素k到標(biāo)記椎骨中心點(diǎn)i的距離.概率密度函數(shù)為p(Ck|Sk).隨機(jī)森林回歸目標(biāo)函數(shù)I(Dj,θ)定義為:

(1)

本文采用連續(xù)形式的微分熵:

(2)

樣本數(shù)據(jù)點(diǎn)D,采用連續(xù)高斯模型密度分布p(Ck),并將b元高斯的微分熵函數(shù)定義為:

(3)

根據(jù)一維輸入與輸出的回增益表達(dá)式推導(dǎo)過(guò)程[14],推廣到多元變量的情況可得:

(4)

其中,Λy為條件協(xié)方差矩陣,可以得到一棵決策樹(shù)的后驗(yàn)概率pt(Ck|Sk),對(duì)T棵樹(shù)的后驗(yàn)概率求平均值,T為決策樹(shù)數(shù)量,并將后驗(yàn)概率平均值作為回歸的結(jié)果.

(5)

因?yàn)闆Q策樹(shù)之間的決策結(jié)果差異較大,隨機(jī)森林中決策樹(shù)的性能越好的決策樹(shù)應(yīng)該占有更高的比重.為解決這個(gè)問(wèn)題,對(duì)隨機(jī)森林中的決策樹(shù)采取加權(quán)運(yùn)算,來(lái)提高準(zhǔn)確率.每棵決策樹(shù)的可信度可由公式(5)計(jì)算得到的后驗(yàn)概率表示.

隨機(jī)森林初始化時(shí),所有樹(shù)的權(quán)值相同,為:

ω0(k)=1/T

(6)

本文提出隨機(jī)回歸森權(quán)重公式為:

(7)

其中,ω(k)表示分類(lèi)森林中第k棵決策樹(shù)的權(quán)重,γ(k)為第k棵決策樹(shù)的回歸中心與專(zhuān)家標(biāo)記點(diǎn)間的距離,γ(k)越大所占權(quán)重越小,通過(guò)迭代優(yōu)化最小絕對(duì)誤差函數(shù)來(lái)訓(xùn)練隨機(jī)回歸森林.

4.1 水平集演化分割

由于水平集分割模型對(duì)初始輪廓敏感且脊柱CT圖像的灰度不均勻分布,可能會(huì)導(dǎo)致陷入局部最小,分割不完全.為了減小誤差,將初始輪廓置于隨機(jī)森林定位的椎骨中心點(diǎn)上并設(shè)定曲線(xiàn)由內(nèi)向外演化.

令I(lǐng)為區(qū)域Ω上的一幅椎骨CT圖像,將分割的初始輪廓L選取在隨機(jī)森林算法確定的椎骨中心點(diǎn)位置,并使用水平集分割算法對(duì)椎骨CT進(jìn)行分割,設(shè)定輪廓由內(nèi)向外膨脹.

將邊緣指示函數(shù)定義為:

(8)

其中Gσ是標(biāo)準(zhǔn)方差為σ的高斯核函數(shù),像素為D(x,y).引用兩個(gè)參數(shù)來(lái)調(diào)整水平集演化的速度及噪聲敏感度,β>0,γ>0,其中β用來(lái)調(diào)整演化速度,γ用來(lái)控制噪聲敏感度.

引進(jìn)參數(shù)后,將新的水平集能量方程定義如公式(9).

E(φ)=αRp(φ)+λLl(φ)+ηAl(φ)

(9)

(10)

使用梯度下降法將公式(7)的能量方程進(jìn)行最小化求解得到水平集演化方程如公式(11):

(11)

通過(guò)逆向有限差分算法,可得:

(12)

根據(jù)公式(10)進(jìn)行水平集函數(shù)φ的演化,直至演化曲線(xiàn)達(dá)到穩(wěn)定狀態(tài)后停止演化.為使分割更為準(zhǔn)確,對(duì)達(dá)到穩(wěn)定狀態(tài)的停止演化的曲線(xiàn)進(jìn)行平滑曲線(xiàn)操作,作為最終的分割輸出結(jié)果,完成椎骨CT圖像分割.

本文方法實(shí)現(xiàn)椎骨CT圖像分割步驟如圖2所示.

圖2 分割步驟Fig.2 Segmentation step

4.2 分割結(jié)果

本文使用公開(kāi)的椎骨CT分割挑戰(zhàn)數(shù)據(jù)集,對(duì)數(shù)據(jù)集中提供的20組脊柱CT圖像共計(jì)1190張CT圖像進(jìn)行分割實(shí)驗(yàn).實(shí)驗(yàn)環(huán)境為Intel Pentium 3.2 GHz CPU、4GHz GPU、MATLAB R2010a.實(shí)驗(yàn)結(jié)果如圖3所示.

在圖3 (a)為初始輪廓椎骨CT的分割結(jié)果;圖3(b)為椎骨中心點(diǎn)位置,其中黑色點(diǎn)為最終確定椎骨中心點(diǎn),周?chē)鷾\色區(qū)域?yàn)榇ㄖ行狞c(diǎn);圖3(c)輪廓分割初始輪廓選取結(jié)果,其中白色方框?yàn)槌跏驾喞?圖3(d)本文方法椎骨分割結(jié)果,深色區(qū)域?yàn)楸疚姆椒ǚ指罱Y(jié)果,白色區(qū)域?yàn)閷?zhuān)家手動(dòng)分割真實(shí)值;圖3(e)為整個(gè)胸椎和腰椎的分割結(jié)果圖.

本文方法對(duì)椎骨可以準(zhǔn)確、穩(wěn)定地分割,相比其他算法對(duì)于任意位置初始輪廓的分割結(jié)果,準(zhǔn)確性具有很大的提高.為了能夠更直觀(guān)的觀(guān)察椎體的分割效果,對(duì)腰椎椎骨L3、L4、L5的分割結(jié)果進(jìn)行了三維重建,如圖4所示.

根據(jù)骰子系數(shù)(DC),對(duì)本文方法進(jìn)行評(píng)估和對(duì)比,其定義如下:

(13)

圖3 椎骨切片分割結(jié)果Fig.3 Vertebral slice segmentation result

其中,Sr是參考體積,Ss是分割體積.將脊柱分為三段(T1-T6,T7-T12和L1-L5)對(duì)分割的平均骰子系數(shù)進(jìn)行計(jì)算,本文方法分割結(jié)果骰子系數(shù)平均為0.936,具體結(jié)果如表1所示,隨機(jī)抽取4組脊柱CT進(jìn)行了分割結(jié)果的骰子系數(shù)統(tǒng)計(jì),結(jié)果如圖5所示.

圖4 腰椎分割結(jié)果三維可視化Fig.4 3D visualization results of lumbar segmentation

T1-T6T7-T12L1-L5All骰子系數(shù)0.9250.9320.9510.936

表2對(duì)本文方法的椎骨分割結(jié)果與其他3個(gè)算法進(jìn)行了詳細(xì)比較.方法1只分割了腰椎區(qū)域,不能分割胸椎區(qū)域;方法1和方法2需要對(duì)椎骨進(jìn)行手動(dòng)定位;方法3和本文方法能夠?qū)ψ倒亲詣?dòng)定位,本文方法分割更準(zhǔn)確.本文方法分割骰子系數(shù)平均為0.936,能夠自動(dòng)定位中心位置并有較高的分割準(zhǔn)確率.

表2 椎骨CT分割結(jié)果對(duì)比
Table 2 Comparison of CT segmentation results of vertebrae

方法椎骨定位分割方法骰子系數(shù)DC方法1[15]手動(dòng)定位統(tǒng)計(jì)形狀模型0.930方法2[16]手動(dòng)定位平均形狀模型0.934方法3[17]自動(dòng)定位平均形狀模型0.930本文方法 自動(dòng)定位WRF-CV 模型0.936

圖5 測(cè)試CT集分割結(jié)果Fig.5 Test CT set segmentation results

5 結(jié)束語(yǔ)

本文提出了基于加權(quán)隨機(jī)森林和CV模型的椎骨CT圖像分割方法.通過(guò)對(duì)不同數(shù)據(jù)進(jìn)行實(shí)驗(yàn)得出,本文提出的算法能夠成功定位椎骨中心位置,能夠有效解決因初始輪廓不正確而造成的過(guò)分割和分割不完全的問(wèn)題,在自動(dòng)定位和準(zhǔn)確性方面具有一定的優(yōu)勢(shì).在之后的研究學(xué)習(xí)中,還需要研究更為通用的椎骨CT圖像分割方法,提高自適應(yīng)性.