陶永鵬,景 雨,頊 聰

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

1 概 述

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

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

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

2 本文方法

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

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

3 SIFT特征

SIFT特征(Scale-invariant feature transform)是一種圖像局部特征提取算法,能夠保持較好的穩(wěn)定性.在描述子的生成過程中,對于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)過程[14],推廣到多元變量的情況可得:

(4)

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

(5)

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

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

ω0(k)=1/T

(6)

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

(7)

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

4.1 水平集演化分割

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

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

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

(8)

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

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

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

(9)

(10)

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

(11)

通過逆向有限差分算法,可得:

(12)

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

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

圖2 分割步驟Fig.2 Segmentation step

4.2 分割結(jié)果

本文使用公開的椎骨CT分割挑戰(zhàn)數(shù)據(jù)集,對數(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),周圍淺色區(qū)域?yàn)榇ㄖ行狞c(diǎn);圖3(c)輪廓分割初始輪廓選取結(jié)果,其中白色方框?yàn)槌跏驾喞?圖3(d)本文方法椎骨分割結(jié)果,深色區(qū)域?yàn)楸疚姆椒ǚ指罱Y(jié)果,白色區(qū)域?yàn)閷＜沂謩臃指钫鎸?shí)值;圖3(e)為整個(gè)胸椎和腰椎的分割結(jié)果圖.

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

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

(13)

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

其中,Sr是參考體積,Ss是分割體積.將脊柱分為三段(T1-T6,T7-T12和L1-L5)對分割的平均骰子系數(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對本文方法的椎骨分割結(jié)果與其他3個(gè)算法進(jìn)行了詳細(xì)比較.方法1只分割了腰椎區(qū)域,不能分割胸椎區(qū)域;方法1和方法2需要對椎骨進(jìn)行手動定位;方法3和本文方法能夠?qū)ψ倒亲詣佣ㄎ?本文方法分割更準(zhǔn)確.本文方法分割骰子系數(shù)平均為0.936,能夠自動定位中心位置并有較高的分割準(zhǔn)確率.

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

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

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

5 結(jié)束語

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