基于水平集和最大穩定極值區域的頸椎椎體分割方法

蘭添才,陳　俊,張怡晨,李翠華*

(1.廈門大學信息科學與技術學院,福建 廈門 361005;2.龍巖學院信息工程學院,福建 龍巖 364000;3.龍巖市第二醫院康復科,福建 龍巖 364000)

近年來,由于長時間面對電腦或手機,頸椎病的患病率逐年增高,發病年齡有明顯年輕化的趨勢[1].X射線、CT及磁共振成像(magnetic resonance imaging,MRI)等影像學檢查已經成為頸椎病診斷的常規輔助手段,相較于CT和MRI,X射線檢查因其最簡便、直觀、經濟且有效,目前仍然是頸椎病早期診斷的主要檢查方法[2].頸椎圖像的精確分割可以輔助醫生判斷病情，量化分析病灶區域,為頸椎病診斷提供可靠的依據,在頸椎病的臨床診斷中有著非常重要的研究價值.然而,頸部是人體幾何特性和運動特性最為復雜的部分之一,由于患者個體的差異性、頸椎結構本身的復雜性和椎體周圍組成物質的非均勻性等,導致難以實現頸椎圖像的精確分割.對頸椎的分割,已有一些學者做了相關研究[3-5].如,Chen等[3]對頸椎形態進行了定量分析,但其采用人工定點分割,依賴于操作者的主觀經驗,且耗時長、重復性低;藍智聰等[4]利用拋物線分別對椎體的四條邊緣進行擬合,進而提取頸椎椎體,其提取過程采用人機交互;趙曉光等[5]在手工截取頸椎區域的基礎上,提出基于邊緣和角點的算法提取椎體左右邊緣,并分別用直線和曲線擬合出椎體的上下邊緣,但其魯棒性有待提高,僅對頸椎的總體排列和其前方的咽后壁呈“直線”型的頸椎有較好的效果.

水平集方法在處理醫學圖像分割問題時具有良好的性能,很多學者在這方面取得了一定的研究成果[6-10].如,梁禮明等[6]將水平集模型應用在眼底圖像血管的分割中取得了良好的效果;Jeffee等[7]和Lu等[8]把水平集方法應用于超聲圖像的子宮頸以及宮頸細胞的分割中,成功提取宮頸以及成功分割出單個宮頸細胞;Wang等[9]提出了一種基于概率圖譜與形狀強度水平集相結合的分割方法并用于肝臟CT圖像的自動分割,取得了類似于手工分割的肝臟邊界;Cheng等[10]提出了一種基于B-Snake模型的血管CT圖像自動分割方法,算法可以分割出類似于手工獲得的血管邊界.醫學圖像分割通常采用多種分割方法的融合[11].水平集方法對于對比度低或者邊緣比較模糊的圖像分割不準確,而頸椎椎體后緣相對于前緣,明顯具有對比度低、邊緣模糊的特點.考慮到后緣具有亮度高、不連續等穩定的局部特征.而對于圖像的局部不變特征提取,最大穩定極值區域(maximally stable extremal regions,MSER)在大部分情況下性能最佳[12].本研究在水平集分割的基礎上,提出采用改進的MSER方法進一步提取頸椎椎體的后緣,從而提取出完整的頸椎椎體.

圖2　自動提取頸椎區域Fig.2 Automatic extraction of cervical regions

1　算法流程及原理

針對頸椎椎體形態規律一致,后緣存在局部穩定特征的先驗,文中提出一種融合水平集和MSER的頸椎椎體分割方法.算法的分割流程如圖1所示.

圖1　分割流程圖Fig.1 The flow chart of segmentation

1.1　頸椎區域的自動提取與圖像增強

頸椎側位X射線圖像來源于拍攝的頭顱定位側位片,如圖2(a)所示.原圖包含整個頭顱的內容,而頸椎病診斷臨床中所關注的區域是七小節頸椎椎體所在的區域.頸椎區域的提取可以避免對全圖進行計算,提高處理速度.如何快速而準確地將頸椎區域從原圖中分離出來,是整個臨床和研究分析過程的首要步驟.趙曉光等[5]采用手工提取頸椎區域費時費力.根據頸椎側位X射線圖像中椎體位置的布局,本文中提出一種基于圖像密集度分布的自動提取算法,根據閾值分割結果的密集度分布來自動提取頸椎區域.具體的提取步驟如下:

1) 降維處理.將圖像統一從1 792像素×2 392像素降維到600像素×800 像素,降低原圖分辨率.

2) 采用最大類間方差法(OTSU)對降維后的圖像進行閾值分割(粗分割),分割結果如圖2(b)所示.

3) 對閾值分割后的結果分別進行水平和垂直方向的投影，投影直方圖如圖2(c)所示.

4) 根據水平和垂直方向的投影直方圖自動獲取閾值,并將圖像大小歸一化為192像素×512像素，頸椎區域自動提取效果如圖2(d)所示.

由圖2(b)可知,對閾值分割后的結果進行垂直方向投影,頸椎部分投影值最高;對分割圖進行水平方向投影,頸椎部分投影具有中間平穩的特性.圖2(c)進一步驗證了這個特征,水平方向投影的結果得到頸椎區域的水平方向的位置,垂直方向的投影結果得到頸椎區域的垂直方向位置.從圖2(d)可知,提取結果完整包括七小節頸椎椎體,且縮小了范圍,排除了上方顱骨和左上方下頜骨等較多的干擾.

受到采集系統和采集條件如光照等諸多因素的影響,采集到的X射線圖像的質量不一致,噪聲較多,這會對后續的圖像處理產生不利的影響,因此有必要對粗分割后的頸椎區域進行圖像增強.文中的圖像增強主要包括以下2種方法:

1) 根據文獻[13],對頸椎區域圖像進行空域上的同態濾波以減少頸椎X射線圖像拍攝時不均勻光照產生的影響.

2) 采用中值濾波對圖像進行去噪，濾波算子為3×1.由先驗知識可知,X射線圖像中多是椒鹽噪聲,而且椎體的邊緣信息特別是椎體后緣的高亮信息是本文中提取頸椎椎體的重要線索.中值濾波的核心思想是圖像中像素灰度值設置為該點某領域窗口內所有像素灰度值的中值,該方法對消除頸椎圖像的椒鹽噪聲非常有效,而且可以較好地保留圖像的邊緣細節,不破壞頸椎椎體后緣的局部高亮等穩定特征.

1.2　基于水平集的頸椎前緣輪廓提取

Chan和Vese[14]提出的C-V模型是一種典型的基于全局區域的水平集分割模型,C-V模型把目標與背景考慮為全域二聚類問題,其缺點主要是不能成功地分割亮度不均勻的圖像.Caselles等[15]提出的測地線活動輪廓(GAC)模型是基于局部邊界的水平集分割模型,缺點是不能成功地分割對比度低及邊緣模糊的圖像,容易造成邊界泄露,且其存在對凹陷處難分割的問題.文獻[16]表明,若用GAC模型中邊界能量項代替C-V模型能量泛函中閉合邊界曲線C的全弧長項,可獲得更好的邊界信息.假設定義域為Ω的圖像I(x,y)被閉合邊界曲線C劃分為內部目標區inside(C)和外部背景區outside(C).融合區域和邊界特征的能量項后,新的能量泛函表示如下[16]:

(1)

其中,ci(i=1,2)表示曲線內部和外部的平均灰度值,λi≥0為能量項權重系數,系數μ>0,g為邊界停止函數.

只有當輪廓線C演化到目標邊界時,能量泛函才能達到最小值.優化式(1),即可得到未知數c1、c2以及最終分割輪廓線C的位置.

分析圖2(d)可知,經過頸椎區域的提取,頸椎椎體形態較相似,每個椎體的前緣清晰完整,亮度分布均勻,椎體之間有間隙,呈現凹凸狀.因此本文中采用C-V模型和GAC模型相結合的方法提取頸椎的前緣輪廓.

1.3　基于MSER的頸椎后緣曲線提取

頸椎椎體后緣X射線圖像實際為椎體后緣與兩側橫突骨皮質的重疊影像,X射線成影后該區域的亮度值較大,形成局部高亮,同時團聚成具有一定面積大小、分布不連續的長條型區域.文獻[12,17]的研究均表明MSER算法在穩定性、仿射不變性和對光照的適應性等方面都優于其他的區域特征提取算法.本文中從頸椎椎體后緣的局部高亮特征出發,利用MSER算法檢測出高亮區域,從而實現頸椎椎體的后緣曲線提取.

1.3.1頸椎MSER區域提取

MSER算法由Matas等[18]于2002年提出,MSER數學定義為:

定義區域D的連通子集Q:對于?p,q∈Q,都存在路徑p,a1,a2,…,an,q,使得pAa1,a1Aa2,…,anAq相鄰,這里ai∈Q,i=1,2,…,n;區域Q的邊界?Q:?Q={q∈DQ:?p∈Q:qAp},即?Q當中的像素不屬于子集Q,但與Q內至少一個像素相鄰.則對于極值區域Q?D,?p∈Q,?q∈?Q,有:

(2)

假設Q1,…，Qi-1，Qi,…為嵌套極值區域的序列,即Qi-1?Qi.若穩定方程:

q(i)=|Qi+△-Qi-△|/|Qi|，

(3)

其中,Qi表示灰度閾值為i時的連通區域,△為灰度閾值的微小增量,q(i)表示區域Qi的面積變化率.若i*處存在局部最小值,則極值區域Qi*即為MSER.

MSER提取過程具體包括以下步驟:1) 對給定的圖像采用桶排序(bucket sort),按照亮度大小對所有像素點進行排序得到灰度值遞增的像素點序列；2) 根據式(2)檢測極值區域；3) 根據式(3)提取MSER；4) 由先驗可知,頸椎椎體后緣多為長條形高亮區域,所以采用矩形擬合方法把問題轉化為對MSER最小外接矩形的分析,再根據外接矩形的參數判斷是否是頸椎椎體后緣區域.

1.3.2噪聲區域分析與剔除

由于整幅圖像中頸椎椎體后緣區域的像素點的像素值有著較大相似性,呈現高亮且不連續的短邊緣特征,所以可認為頸椎的每一節椎體的后緣高亮區域可能對應一個MSER,因為不連續也可能對應多個MSER.但是因為頸椎背景復雜,經過粗分割后的頸椎區域還包括上方的部分顱骨、左上方部分下頜骨以及后緣周圍的棘突、橫突等,這些區域也可能存在高亮部分,在對圖像進行MSER提取后,除了在椎體后緣所在區域提取到MSER外,在上述這些非目標區域可能也會提取到一些MSER.為了有效剔除這些噪聲區域的MSER,同時又能最大限度地保留椎體后緣上的MSER,文中提出以下策略:

策略一:剔除不符合高度位置的MSER.從椎體結構的先驗可知,第1,2節椎體結構特殊,故本文中主要提取第3～7節頸椎椎體,所以特征點也主要取自第3～7節頸椎椎體,第2節椎體以上的MSER作為噪聲區域剔除.雖然由于個人體質、形變和拍攝角度導致椎體位置不盡相同,但是各節椎體高度大致相同,共7節椎體,本文中取整體椎體高度的2/7作為閾值剔除第1，2節椎體中的MSER.

策略二:剔除不滿足橫向間距的MSER.由頸椎結構先驗可知,頸椎椎體后緣的縱向偏移應在±15°范圍內,所以相鄰椎體之間的橫向距離也應在較小距離范圍內.如圖3(a)所示,假設(x1,y1)、(x2,y2)為相鄰椎體MSER矩形的中心點坐標,則兩個中心點的橫向距離|x2-x1|應在一個較小的閾值范圍內.按縱軸方向,通過計算相鄰MSER對應矩形的中心位置的橫向偏移,可剔除不滿足條件的MSER.具體操作如下:首先掃描圖像中經策略一處理之后的所有經過區域擬合后的矩形區域,分別計算每個矩形的中心點位置坐標,然后按照中心點位置的縱坐標大小對矩形進行排序.假設共有k個矩形Rk,其對應的中心點位置坐標為(xi,yi),i=1,2,…,k，則相鄰中心點的平均距離為

(4)

圖3　噪聲區域剔除策略示意圖Fig.3 Schematic diagram of the noise regions excluding strategy

以d0再加上一個增量ζ0作為閾值,用以剔除不符合條件的噪聲MSER.根據大量實驗,本文中設定ζ0為0.2 mm.

(5)

5) 結合椎體的前后緣寬度先驗知識,以d1加上一個增量ζ1作為閾值,剔除不符合條件的MSER區域.根據大量實驗,本文中設定ζ1為0.7 mm.

1.3.3頸椎后緣曲線提取

通過提取椎體的MSER區域并排除噪聲區域后,提取所有目標區域MSER的矩形中心點,采用最小二乘法擬合出椎體后緣曲線,但由于排除噪聲區域后的MSER主要集中在頸椎第3～7節椎體,而且MSER 特征形狀不規則,可能會導致后緣曲線的提取存在誤差.為此,本文中采用對后緣曲線進行二次擬合的方法,首先在擬合出的后緣曲線上獲取曲線的頂點位置坐標,并參照臨床中常用的Boden氏測量方法[19]提取第2節頸椎椎體后上緣和第7節頸椎后下緣位置上的兩點坐標,根據這兩點坐標對曲線進行二次擬合,得到更準確的后緣曲線.

1.4　算法流程

本文中頸椎椎體分割算法(簡稱本文算法)的具體流程如下:

1) 采用基于圖像密集度分布的分割算法自動提取頸椎區域;

2) 對提取的頸椎區域進行圖像增強;

3) 采用C-V和GAC相結合的水平集算法(CV-GAC)提取頸椎椎體前緣輪廓;

4) 采用MSER算法提取頸椎區域的MSER;

5) 對提取的前緣輪廓進行跟蹤并細化處理;

6) 結合椎體結構剛性特征,通過文中提出的3種策略剔除噪聲區域的MSER;

7) 采用最小二乘法對剔除噪聲區域后的所有MSER進行擬合,得到椎體后緣曲線;

8) 參照Boden氏測量方法,分別提取第2節頸椎椎體后上緣(椎體后緣中偏上)和第7節頸椎后下緣位置上的兩點坐標;

9) 結合人工提取的兩點坐標,對后緣曲線采用最小二乘法進行二次擬合;

10) 融合前緣輪廓和后緣曲線,提取頸椎椎體區域.

2　實驗結果與分析

2.1　實驗環境與實驗數據

本文中的圖像分割實驗是在Windows系統下,利用MATLAB R2014b實現的.Lenovo Thinkpad 筆記本具體配置為 CPU Intel Core(TM) i5 2.60 GHz,4 GB RAM,并采用 SPSS 22.0 軟件進行統計學分析.

實驗數據來源:本文中實驗的所有頸椎X射線圖像均由龍巖市第二醫院放射科同一組專業技師完成,攝片條件相同(機器所設置的參數),且拍攝頭顱側位圖像時要求受試者頸部姿勢統一.實驗選取頸椎病臨床診斷中常見的生理曲度正常、反弓、變直以及旋轉等4種類型X射線頸椎側位圖像共計170例進行算法驗證.

2.2　實驗結果分析

根據本文算法流程,首先要提取椎體的前緣輪廓,圖4分別給出了頸椎生理曲度正常、反弓、變直以及旋轉4種常見類型的X射線原圖(圖4(a)),CV-GAC提取的頸椎前緣輪廓效果圖(圖4(b)),專家組醫生手工分割的頸椎前緣效果圖(圖4(c)).實驗中基于區域的能量函數權重取常用值,即λ1=λ2=1;在實驗中發現μ的取值對椎體右側輪廓影響較大,對左側輪廓影響較小,因椎體右側結構導致右側背景較復雜.而本研究關注的是椎體的左側輪廓,實驗中μ的值取為500.

從4(b)可知,本文中提出的自動提取頸椎區域的方法簡單有效,提取的頸椎區域完整;采用CV-GAC方法提取的頸椎椎體前緣輪廓曲線非常貼近頸椎真實前緣,幾乎與頸椎椎體的前緣重合.與圖4(c)比較可知,采用CV-GAC方法提取的頸椎椎體前緣接近專家組手工分割的結果.對于存在較嚴重病變,如旋轉扭曲的頸椎病變圖像,雖然有誤檢,但實驗中發現這種誤檢多發生在旋轉圖像的第2節椎體,是因為旋轉使舌骨在該處形成重影所致,而其他椎體由于具有相對穩定的剛性結構特征,同樣能夠得到貼近椎體前緣的較完整輪廓曲線,這將為后續頸椎椎體后緣曲線的準確提取起關鍵作用.

圖4　CV-GAC方法與手工分割方法提取頸椎前緣輪廓Fig.4 Front edge of cervical contour extracted by CV-GAC and manual segmentation method

其次,要提取頸椎區域的MSER并排除噪聲區域的MSER,這是后緣曲線提取的關鍵.圖5(a)是4種類型的頸椎區域MSER提取效果圖；圖5(b)所示為剔除噪聲區域MSER后的效果圖,圖中被矩形框住的區域表示選中區域,沒有被矩形框住的表示被剔除的噪聲區域MSER.

從圖5(a)可知,各類頸椎圖像中第2～7節椎體后緣目標區域的絕大部分高亮MSER均能被有效檢出,盡管旋轉圖像較特殊,仍能有效檢出.但同時有較多噪聲區域的MSER被檢測到,且噪聲區域的MSER主要集中在圖像上部顱骨等處,椎體周圍檢測出的噪聲區域的MSER較少.對于圖像上部顱骨等處的噪聲區域的MSER采用策略一容易剔除;而椎體周圍檢測出的噪聲區域的MSER雖然較少,但它對后續的椎體后緣曲線擬合影響較大,必須剔除.本文中剔除噪聲區域MSER的3種策略都是基于形狀結構信息而提出.實驗表明,采用本文中提出的3種策略能有效剔除噪聲區域的MSER,同時又能完整保留椎體后緣高亮區域的MSER.剔除效果如圖5(b)所示.

圖5　頸椎區域的MSER提取與噪聲區域的剔除Fig.5 MSER extraction of cervical reguons and exclusion of the noise regions

圖6　本文算法與手工分割結果比較Fig.6 Results comparison of proposed segmentation and the manual segmentation

圖6(a)是采用最小二乘法對排除噪聲后剩下的所有MSER區域進行擬合,得到后緣曲線,并與前緣輪廓進行融合后的結果.為使椎體后緣擬合得更準確,便于臨床診斷數據的提取,參照臨床中常用的Boden氏測量方法,分別提取第2節椎體后緣中偏上位置和第7節椎體后下緣位置上的兩點坐標,并對椎體后緣曲線進行二次擬合,擬合的效果如圖6(b)所示.圖6(c)給出了專家組醫生采用手工分割方法提取的頸椎后緣曲線效果圖.根據文獻[19]中的診斷標準,從圖6(b)可知,正常、反弓和變直3種類型的椎體后緣曲線擬合的效果都很好,與實際的頸椎椎體的后緣都非常的貼近;雖然旋轉類型圖像本身特殊,但從擬合的角度看也比較貼近.與6(c)比較可知,擬合結果接近專家組醫生手工分割的結果.

從表1也可知,本文算法與專家組手工分割提取的曲度間不存在顯著性差異(p>0.05),統計結果表明本文算法提取的后緣曲線接近專家組醫生手工分割提取的結果.

為進一步驗證本文算法的有效性,根據陳銀海等[20]對曲度異常的頸椎分類標準,將頸椎圖像分為正常(7 mm≤曲度值≤17 mm)、反弓(曲度值≤1 mm)和變直(1 mm<曲度值<7 mm)3類.并根據之前提取的曲度值,從170例中篩選出區分度較好的正常圖像55例、反弓圖像17例以及變直圖像28例(共100例),采用SPSS 22.0再次進行統計學分析比較,統計結果如表2所示.從表2可知,采用本文算法提取的3種類型的頸椎后緣曲線的曲度值均與手工分割結果非常的接近,這進一步表明本文算法提取的后緣曲線接近專家組醫生手工提取的結果,特別是對正常、反弓及變直的頸椎圖像可獲得理想的效果.

表1　本文算法和手動分割測得的頸椎曲度

表2　3種類型頸椎不同方法測得的頸椎曲度

注：橫向比較有顯著性差異(p<0.05);縱向比較無顯著性差異(p>0.05).

3　結　論

頸椎側位X射線圖像是目前頸椎診斷中最常用的影像學檢查方法,頸椎椎體的分割在X射線頸椎圖像數據測量中起著關鍵的作用.但是因為測量方法受限,臨床實際操作中普遍采用人工分割為主,人工分割主觀性較強、效率低,對此本文中提出了基于水平集和MSER的頸椎椎體分割方法.從實驗結果看,該方法能完整提取頸椎椎體,對正常、反弓及變直的頸椎椎體特別有效,尤其是后緣曲線的提取與專家手工提取的結果非常接近,可為頸椎曲度的測量、相鄰椎體屈曲度的測量以及頸椎骨齡判定中的測量等實際的測量操作提供較準確且客觀的數據.因此,本文算法在頸椎的輔助診斷中具有一定的臨床實用價值.

參考文獻：

[1]張少群,李義凱.頸椎病研究的歷史沿革[J].中國康復醫學雜志,2016,31(11):1273-1276.

[2]陳超.探討不同的影像學方法診斷頸椎病的臨床價值[J].檢驗醫學與臨床,2013,10(20):2694-2695.

[3]CHEN L L,LIN J X,XU T M,et al.The longitudinal sagittal growth changes of maxilla and mandible according to quantitative cervical vertebral maturation[J].Journal of Huazhong University of Science and Technology (Medical Sciences),2009,29(2):251-256.

[4]藍智聰,陳莉莉,許向陽,等.計算機輔助頸椎分析系統識別頸椎標志點的準確性和重復性探討[J].臨床口腔醫學雜志,2010,26(9):526-530.

[5]趙曉光,林久祥,王晴竹.頸椎影像計算機自動識別骨齡診斷系統的建立[EB/OL].(2013-09-03)[2017-06-17].http:∥www.paper.edu.cn/releasepaper/content/201309-17.

[6]梁禮明,黃朝林,石霏,等.融合形狀先驗的水平集眼底圖像血管分割[J/OL].計算機學報，2016.(2016-11-25).[2017-06-17].http:∥www.cnki.net/kcms/detail/11.1826.TP.20161125.2311.002.html.

[7]JEFFEE A I B,PAHL C,ABDULJABBAR H N,et al.Cervical segmentation in ultrasound image using level-set algorithm[C]∥Advances in Biomedicine and Health Science.Johor:[s.n.],2013:25-30.

[8]LU Z,CARNEIRO G,BRADLEY A P.An improved joint optimization of multiple level set functions for the segmentation of overlapping cervical cells[J].IEEE Transactions on Image Processing,2015,24(4):1261-1272.

[9]WANG J K,CHENG Y Z,GUO C Y,et al.Shape-intensity prior level set combining probabilistic atlas and probability map constrains for automatic liver segmentation from abdominal CT images[J].International Journal for Computer Assisted Radiology and Surgery,2016,11(5):817-826.

[10]CHENG Y Z,HU X,WANG J,et al.Accurate vessel segmentation with constrained B-Snake[J].IEEE Transactions on Image Processing,2015,24(8):2440-2455.

[11]王陽萍,杜曉剛,趙庶旭,等.醫學影像圖像處理[M].北京:清華大學出版社,2012:69-109.

[12]MIKOLAJCZYK K,TUYTELAARS T,SCHMID C,et al.A comparison of affine region detectors[J].International Journal of Computer Vision,2005,65(1/2):43-72.

[13]聞莎,游志勝.性能優化的同態濾波空域算法[J].計算機應用研究,2000,17(3):62-65.

[14]CHAN T F,VESE L.Active contours without edges[J].IEEE Transactions on Image Processing,2001,10(2):266-277.

[15]CASELLES V,KIMMEL R,SAPIRO G.Geodesic active contours[J].International Journal of Computer Vision,1997,22(1):61-79.

[16]TAO W B,TAI X C.Multiple piecewise constant with geodesic active contours(MPC-GAC) framework for interactive image segmentation using graph cuts optimization[J].Image and Vision Computing,2011,29(8):499-508.

[17]MIKOLAJCZYK K,TUYTELAARS T,SCHMID C,et al.A comparision of affineregion detectors[J].International Journal of Computer Vision,2005,65(1):43-72.

[18]MATAS J,CHUM O,URBAN M,et al.Robust wide baseline stereo from maximally stable extremal regions [C]∥Proceedings of the British Machine Vision Conference.Cardiff:[s.n.],2002:384-393.

[19]王濤,周理乾,孫孟錕,等.6種頸椎曲度測量方法的可信度及可重復性比較[J].中國脊柱脊髓雜志,2015,25(4):323-327.

[20]陳銀海,姚紅華,楊忠.頸椎曲度的X線測量在頸椎病康復評定中的應用價值[J].中國康復,2007,22(3):156-158.