王 浩，邵 堃，霍 星，楊 鵬，檀結慶

?

一種基于極大特征點的三維椎骨分割方法

王 浩1，邵 堃1，霍 星2，楊 鵬1，檀結慶2

(1. 合肥工業(yè)大學計算機與信息學院，安徽 合肥 230000；2.合肥工業(yè)大學數(shù)學學院，安徽 合肥 230000)

椎骨的精確分割對于椎骨形態(tài)學研究和脊柱疾病的診斷和治療有重要意義。通過對正常人脊柱CT圖序列的變化規(guī)律進行研究，提出了一種基于CT圖像序列并利用椎骨面積變化規(guī)律進行分割的椎骨分割算法。該方法通過對預先處理后的CT圖像序列進行椎骨區(qū)域面積統(tǒng)計，找出用于分割的顯著極大特征點，并利用連續(xù)圖像相似性篩選出椎骨實際分割點，最后從序列中提取圖像并進行三維重建。實驗表明，該算法對正常人體腰椎和胸椎下部的椎骨CT圖像序列有良好地分割效果，自動化程度較高。對脊柱形態(tài)學研究和矯正手術模擬有重要意義。

CT圖像序列；椎骨分割；面積變化規(guī)律；三維重建

隨著計算機圖像技術的不斷發(fā)展，當今社會數(shù)字圖像處理技術在各個領域都有重要地應用，展現(xiàn)在醫(yī)學領域的應用也越來越廣泛。醫(yī)學影像設備的不斷升級，三維造型技術的不斷進步，數(shù)字診療技術的不斷提升，都離不開醫(yī)學圖像處理技術的發(fā)展。如今，醫(yī)學圖像處理已成為計算機科學和生物醫(yī)學工程領域研究的熱點。

椎骨分割一直是醫(yī)學圖像處理領域中研究的重點和熱點。一個良好的椎骨分割結果對人們研究椎骨形態(tài)、計算機輔助診斷及治療等都具有重大的意義[1]。但是，由于脊柱空間結構較為復雜、各椎骨相互勾連，存在較大的局部噪聲[2]，使得椎骨分割成為一項極具挑戰(zhàn)性的任務，同時也吸引了大批研究者從事相關問題的研究。

近年來，隨著三維手術導航技術的開發(fā)與應用，如何從整體三維脊柱模型中快速分割出完整的單塊脊椎骨并生成三維模型也成為研究熱點之一。

當今醫(yī)學圖像分割的方法主要有：基于圖像、紋理和邊緣3種[3]。目前椎骨分割方法大致有：基于閾值和區(qū)域增長的分割方法[4]；基于主動輪廓模型檢測椎骨區(qū)域的分割方法[5]；基于圖割與均值漂移算法，利用椎骨中心線進行分割的方法[6]；基于深度學習的椎骨定位與分割方法[7]；基于統(tǒng)計模型和剛性配準的分割方法[8]；基于標準脊柱邊緣輪廓檢測框架的分割方法[9]等?？蓪⑸鲜龇指罘椒w為2類：①從CT或MR三維圖像的矢狀面入手，主要集中于椎體的定位、識別與分割；②從主要軸心面入手，在二維圖像上將椎骨區(qū)域分割出來。其中，第一類方法得到的分割結果往往不夠精確，椎骨模型的完整性較差；第二類方法雖然在二維圖像上有較好的區(qū)域分割結果，但往往沒有將椎骨進行分離。

由此可見，目前椎骨分割方法，仍存在以下不足：①所用的算法和數(shù)據(jù)過于復雜，導致時間復雜度很高，效率較低；②實現(xiàn)分割的限制條件較多，例如只能分割腰椎區(qū)域的椎骨；③自動化程度較低，需要的人工干預較多。因此需要對分割方法進一步進行優(yōu)化[3]。

通過研究人體脊柱CT圖序列，發(fā)現(xiàn)在正常脊柱圖像序列中，椎骨區(qū)域面積的大小呈現(xiàn)周期性變化的規(guī)律?；诖吮疚奶岢隽艘环N從正常人體脊柱CT圖序列中快速分割完整椎骨圖像的方法，利用正常脊柱的CT斷層掃描圖，經(jīng)過預處理去除大部分干擾因素后得到二值圖像；然后對預處理后圖像序列進行面積統(tǒng)計和數(shù)據(jù)處理，并根據(jù)面積變化曲線尋找其波谷特征點，利用前后相鄰圖片的面積相似度檢測及篩選這些特征點；最后得到椎骨分割點并利用分割點提取圖像數(shù)據(jù)進行三維重建，得到單塊脊椎骨模型。

1 分割算法

1.1 流程、算法概述

本文的分割方法分為預處理、椎骨區(qū)域面積檢測、波谷特征點提取與篩選和單塊椎骨圖像提取。4個步驟為：①將CT圖像轉(zhuǎn)換成bmp圖像并進行基礎處理，凸顯椎骨部分的圖像特征；②計算圖像序列中每張圖像的椎骨區(qū)域的面積，生成面積變化曲線，提取曲線上的波谷特征點并進行篩選，確定最終的椎骨分界點；③利用分界點在圖像序列中進行圖像提取，得到單塊椎骨的圖像序列；④進行三維建模，得到單塊椎骨模型。

1.2 圖像預處理

由于CT圖像不易直接進行處理，本文對DICOM醫(yī)學圖像的格式和操作進行了調(diào)研與分析[10-13]并確定了處理方法。首先將CT圖轉(zhuǎn)換為容易處理的bmp圖像，同時將一些參數(shù)以文本的形式保存下來。

然后對bmp灰度圖像序列進行二值化，以增強骨組織區(qū)域特征，為后續(xù)分割提供基礎，并對圖像進行中值濾波，其選取方法為

其中，為以當前像素為中心的滑動窗口；和均為中的相對像素坐標；和為圖像中實際的像素坐標；為映射到滑動窗口的像素值；為圖像的實際像素值。濾波原理是先將滑動窗口覆蓋的圖像中的像素按照大小進行排序并得到中值，然后讓窗口中心位置的像素等于計算出的中值。中值濾波對椒鹽噪聲的抑制效果好，在抑制隨機噪聲的同時能有效保護邊緣，減少模糊。中值濾波窗口大小應為奇數(shù)，本文選用9×9的窗口。圖1為圖像進行過預處理的結果。

圖1 經(jīng)過中值濾波后得到的圖像

為了去除大部分由肋骨帶來的干擾，中值濾波后，使用矩形框?qū)⒆倒撬趨^(qū)域框選出來，完成預處理操作。

1.3 極大特征點提取

經(jīng)過預處理可以得到一系列包含椎骨主體的二值圖像，并對其進行觀察，研究發(fā)現(xiàn)正常的人體脊柱表現(xiàn)的變化規(guī)律為：圖像序列中的椎骨骨質(zhì)區(qū)域的面積呈現(xiàn)波浪形變化，并且當圖像處于兩塊椎骨過渡的位置時，其中骨質(zhì)區(qū)域的面積在相鄰幾張圖片中是最小的。利用此規(guī)律，可以在面積變化曲線上尋找波谷點作為極大特征點，其包含了所有的椎骨分割點；接下來運用相鄰圖像間的相似性，從特征點中篩選出2塊椎骨的分割點；最后，以該分割點為界，從整體圖像序列中提取出單塊椎骨的圖像，并重建成三維模型。

1.3.1 骨質(zhì)區(qū)域面積計算

在含有椎骨的二值圖像中，白色區(qū)域即為骨質(zhì)區(qū)域，但由于骨質(zhì)區(qū)域的不規(guī)則形狀，無法直接計算其面積，但可以通過統(tǒng)計像素灰度值為255的像素點個數(shù)作為骨質(zhì)區(qū)域的面積值。本文以此計算圖像中骨質(zhì)區(qū)域面積的大小。

通過遍歷序列中的每一幅圖像，計算圖像中白色像素點的個數(shù)，最終得到所有圖像的骨質(zhì)區(qū)域的面積值。圖2是正常脊柱CT圖序列得到的脊柱面積變化曲線。

1.3.2 計算曲線波谷特征點

經(jīng)過對椎骨面積變化曲線的觀察和驗證可知：所有的椎骨分割點都處在椎骨面積變化曲線的波谷點上。曲線上的波谷特征點可表示為

其中，為曲線上的橫軸坐標；()為面積值，波谷特征點即面積值小于其前后面積值的那些點。

由于存在噪聲，定位不準，對尋找波谷特征點有一定干擾，為此設置了一個步長間隔，代表2個波谷點之間最小的距離。初步過濾掉一些干擾點，加快計算的速度。圖3是識別波谷點的結果。

圖3 椎骨面積變化曲線中尋找波谷點

1.3.3 篩選波谷特征點

由于并不是所有波谷點都是椎骨分割點，所以需要對波谷特征點進行篩選，為此采用了一種基于面積相似度的方法[14]。該方法通過比較每個波谷點代表的圖像b和其前后相鄰圖像a和c中的椎骨面積，來檢測波谷特征點是否為最終的椎骨分割點。步驟如下：

(1) 根據(jù)圖像序列的間隔距離大小，確定需要比較的3張圖像之間的圖像間隔數(shù)量。若圖像間隔距離大于1，則取波谷點前后相鄰的3張圖像進行比較；若間隔小于1，則取波谷點前后相隔一幅圖像的3張圖像進行比較，其原因是若圖像之間的間距太小，則前后相鄰圖像之間的相似度很高，不利于篩選。然后，分別計算圖像a和c中椎骨區(qū)域area_a與area_c和圖像b的椎骨區(qū)域area_b的重合區(qū)域的大小，得到重合區(qū)域的大小分別占圖像a和c的椎骨區(qū)域的比例，即

(2) 得到2個比例后，利用閾值來判斷1和2是否符合條件，低于說明此處圖片相似度低，處于2塊椎骨的過渡位置，則需要保留此點；反之說明相似度高，還在一塊椎骨內(nèi)部，需去除此點。這里閾值是根據(jù)圖像之間的間距來決定的，經(jīng)測試發(fā)現(xiàn)，2幅圖像間隔越小，值越大才能得到良好的分割效果，根據(jù)多組圖像的分割測試結果，擬合了取值的公式，其中，為序列中圖像之間的間隔距離，即

圖4是最終得到的分割點的示意圖。

圖4 波谷點篩選后的結果

1.4 圖像提取和三維模型重建

得到圖像序列中椎骨之間交界處的分割點后，利用分割點抽取分割點之間的椎骨圖像區(qū)域，最后使用Marching-Cubes算法，對分割后提取出的圖像序列進行三維重建，得到分割后的三維椎骨模型。

2 實驗結果

2.1 正常椎骨分割情況

向醫(yī)院申請獲取了幾組DICOM切片數(shù)據(jù)，并選取了相對正常脊柱的一組數(shù)據(jù)，利用本文方法進行分割操作，圖5是分割所得的結果。對分割后的圖像和原始圖像分別進行三維重建，并將單塊椎骨模型重新組合成完整脊柱，得到了分割結果與實際結果的對比。

圖5(a)中每種顏色為一個分割結果，每2條直線之間是一塊椎骨。從建模的結果來看，發(fā)現(xiàn)在胸椎上部T6之前，分割的效果不是很理想，分割點將一塊椎骨分成了上下兩部分或者未找到正確分割點導致兩塊椎骨被看成一塊。但在胸椎下部和腰椎區(qū)域，分割的效果明顯更好，能夠正確找到分割點并將椎骨完整分割出來，其原因如下：由于在胸椎上部區(qū)域，各椎骨之間間距很小，且有肋骨干擾，所以僅靠面積相似度難以準確判斷波谷點是否為椎骨分割點，而胸椎下部和腰椎區(qū)域，各椎骨之間的間距變大，且肋骨干擾逐漸減小，所以面積相似度可以作為判斷是否為分割點的依據(jù)。

圖5 正常脊柱的分割建模效果

2.2 在椎骨數(shù)據(jù)集上的分割情況

B. Glocker等人為本研究提供了一些開放的脊柱CT圖像數(shù)據(jù)集，其圖像的清晰程度較從醫(yī)院獲取的圖像相比稍有遜色，但不影響使用此方法進行分割。圖像效果對比如圖6所示。

圖6 2種數(shù)據(jù)圖像質(zhì)量情況對比

開放數(shù)據(jù)集中包含138位病人的部分脊柱序列圖，從中選取了45組相對正常的脊柱序列圖并利用本文方法對其進行了分割，并按照數(shù)據(jù)中有無腰椎數(shù)據(jù)將45組數(shù)據(jù)分為2組，一組含有腰椎，另一組不含腰椎。并按照每組數(shù)據(jù)中實際含有椎骨數(shù)量和實驗分割出的椎骨數(shù)量對實驗后的數(shù)據(jù)進行統(tǒng)計。圖7和圖8是2組數(shù)據(jù)的統(tǒng)計直方圖。

圖7 含有腰椎的實驗組分割情況

圖8 不含腰椎的實驗組分割情況

圖7和圖8中藍色代表實際的椎骨數(shù)量，白色代表實驗分割出的椎骨數(shù)量。從圖中可以看出，含有腰椎的實驗組分割椎骨數(shù)量與實際椎骨數(shù)量差異不大，而不含腰椎的實驗組差異較大。本文用皮爾遜相關系數(shù)量化了2組實驗組中實際含有的椎骨數(shù)量和分割出的椎骨數(shù)量之間的相關性，含有腰椎及不含腰椎的實驗組相關系數(shù)分別為0.961 9，0.759 0。由此可見，含有腰椎的比不含腰椎的實驗組的分割情況要好得多。

2.3 在側(cè)彎脊柱上的分割結果

為了驗證本方法的適用范圍，本文挑選了脊柱側(cè)彎患者的CT圖序列，此患者在T4-L2部位的脊柱側(cè)彎比較嚴重。利用本文方法對其脊柱CT序列進行了分割，并將分割結果與未分割的相應部分脊柱進行了對比，結果如圖9所示。

圖9 2種數(shù)據(jù)圖像質(zhì)量情況對比

圖9(a)是分割結果，每種顏色代表一個分割結果，圖9(b)是未分割的整體脊柱。由圖可知，本文的方法在彎曲嚴重的脊柱區(qū)域分割效果不理想，沒有分割出一塊完整的椎骨，只有在相對位置正常的L2-L4腰椎部分，才分割出了相對完整的椎骨。其原因是：彎曲脊柱表現(xiàn)在二維CT圖上的形狀復雜多變，且面積變化規(guī)律不明顯，因此只有在位置相對正常的椎骨處才能得到良好的分割效果。

2.4 與同類分割方法的比較

經(jīng)文獻查閱，發(fā)現(xiàn)文獻[3]與本文方法和目的均較相似。文獻[3]第一步首先對包含椎骨的CT圖像進行聚類并進行分割，然后對聚類后的區(qū)域進行遍歷，遍歷時分為圖像內(nèi)遍歷和間遍歷，圖像內(nèi)遍歷得到單張圖像中的椎骨區(qū)域，利用圖像間遍歷來確定椎骨的邊界；第二步將標記后的椎骨區(qū)域重建為3D粗模型，然后利用3D插值，將粗模型細化得到最終的分割效果。

與本方法相比，文獻[3]在第一步對圖像進行處理時，聚類并進行3D圖像椎骨分割的時間復雜度較高，隨后進行3D插值時，算法時間復雜度也較高；而本方法僅對圖像進行了基礎的二值化以及濾波操作，然后利用像素統(tǒng)計的方法進行單塊椎骨圖像抽取，并且重建后的三維數(shù)據(jù)不需要進一步操作，所以本方法在時間效率上有一定優(yōu)勢。

因文獻[3]將圖像進行了聚類，并以聚類區(qū)域為操作單位，而本文方法以整張圖像作為操作單位，因此在椎骨分割的精度、完整性以及方法適用性上，文獻[3]的方法較本文方法更好。

3 結束語

本文提出的是一種基于波谷特征點的椎骨自動分割方法，利用椎骨區(qū)域在CT圖像序列中的變化規(guī)律，用檢測面積變化曲線中波谷特征點的方式對椎骨分割點進行定位。經(jīng)測試發(fā)現(xiàn)，該方法能夠有效地應用于相對正常的脊柱CT切片圖序列并對其中的椎骨進行提取，對于正常脊柱的胸椎下部以及腰椎部分有較好地分割效果，能夠分割出完整的單塊脊椎骨圖像并生成三維模型。相比于傳統(tǒng)的手工分割方法，本文方法大大提高了分割的效率且避免了繁瑣的勞動，能夠良好地分離正常脊柱中的椎骨。分割出的圖像經(jīng)三維重建后的模型表面光滑且數(shù)據(jù)量少，在脊柱分析、手術模擬、有限元分析等方面有著良好地應用價值。

但是本文方法在正常脊柱的胸椎上部以及側(cè)彎脊柱中的彎曲部分上分割效果較差，原因是這些區(qū)域中椎骨區(qū)域形狀復雜，導致前后椎骨面積變化不夠明顯，僅從面積上難以確定椎骨間的分割點。并且本文方法只能在豎直方向上分離椎骨圖像，對于傾斜狀態(tài)的椎骨不能完整地提取。

針對本文方法存在的不足，在后續(xù)的研究中，將擬利用空間鄰域算法，以得到的椎骨分割點為基礎，對3D圖像進行處理，以提高椎骨分割的完整性以及對彎曲脊柱的適用性。

(感謝：B. Glocker等人為本研究提供了一些開放的脊柱CT圖像數(shù)據(jù)集!)

[1] BAZILA, MIR A H. Segmentation of lumbar intervertebral discs from spine MR images [C]//2014 Innovative Applications ofComputational Intelligence on Power, Energy and Controls with Their Impact on Humanity (CIPECH). New York: IEEE Press, 2014: 85-91.

[2] ZHAO L K, ZHENG S Y, WEI H T. Adaptive active contour model driven by global and local intensity fitting energy for image segmentation [J]. Optik - International Journal for Light and Electron Optics, 2017, 140: 908-920.

[3] COURBOT J B, RUST E, Monfrini E. Vertebra segmentation based on two-step refinement [J]. Journal of Computational Surgery, 2016, 4(1): 1.

[4] KIM Y, KIM D. A fully automatic vertebra segmentation method using 3D deformable fences [J]. Computerized Medical Imaging and Graphics: the Official Journal of the Computerized Medical Imaging Society, 2009, 33(5): 343-352.

[5] ATHERTYA J S, SARAVANA K G. Automatic segmentation of vertebral contours from CT images using fuzzy corners [J]. Computers in Biology and Medicine, 2016, 72: 75-89.

[6] CHENG E K, LIU Y X, WIBOWO H, et al. Learning-based spine vertebra localization and segmentation in 3D CT image [C]//IEEE 13th International Symposium onBiomedical Imaging (ISBI 2016). New York: IEEE Press, 2016: 160-163.

[7] SUZANI A, RASOULIAN A, SEITEL A. Deep learning for automatic localization, identification, and segmentation of vertebral bodies in volumetric MR images [EB/OL]. (2015-03-18)[2018-06-29]. https://doi.org/10.1117/12.2081542.

[8] RASOULIAN A, ROHLING R, ABOLMAESUMI P. Lumbar spine segmentation using a statistical multi-vertebrae anatomical Shape+Pose model [J]. IEEE Transactions on Medical Imaging, 2013, 32(10): 1890-1900.

[9] KOREZ R, IBRAGIMOV B, LIKAR B. A Framework for automated spine and vertebrae interpolation-based detection and model-based segmentation [J]. IEEE Transactions on Medical Imaging, 2015, 34(8): 1649-1662.

[10] LI C, MANTELL S, POLLA D. Design and simulation of an implantable medical drug delivery system using microelectro-mechanical systems technology [J]. Sensors and Actuators A Physical, 2001, 94(1): 117-125.

[11] ESCOTT E J, RUBINSTEIN D. Free DICOM image viewing and processing software for your desktop computer: What’s available and what it can do for you [J]. Radiographics, 2003, 23(5): 1341-1357.

[12] ROSSET A, SPADOLA L, RATIB O. OsiriX: An open-source software for navigating in multidimensional DICOM images [J]. Journal of Digital Imaging, 2004, 17(3): 205-216.

[13] CHUNG H COBXAS D, BIRDSELL L, et al. Automated segmentation of muscle and adipose tissue on CT images for human body composition analysis [C]// Medical Imaging 2009: Visualization, Image-Guided Procedures, and Modeling.Bellingham: International Society for Optics and Photonics, 2009: 72610K-72610K-8.

[14] 姚沛延, 趙宏安, 周明全, 等. 基于CT圖像的脊椎骨組織分割[J]. 計算機工程與應用, 2009, 45(27): 237-239.

A Method of 3D Vertebra Segmentation Based on Maximum Feature

WANG Hao1, SHAO Kun1, HUO Xing2, YANG Peng1, TAN Jie-qing2

(1. School of Computer Science and Engineering, Hefei University of Technology, Hefei Anhui 230000, China; 2. School of Mathematics, Hefei University of Technology, Hefei Anhui 230000, China)

Accurate segmentation of vertebrae has become more and more important for the diagnosis and treatment of spinal diseases. After studying the change rule of some spinal CT image sequences of normal people, we proposed a new vertebral segmentation algorithm based on vertebral area change regulation. This method first calculates the area of vertebral region in the image after pre-processing, then finds out the significant features for segmentation, uses image similarity of continuous images to select the actual segmentation points, and finally extracts vertebral images from image sequence for 3D reconstruction. The innovation of the method is that a new segmentation algorithm based on the maximum feature and correlation of CT image sequence is proposed. The experimental results show that the algorithm has a good segmentation effect on the CT image sequence of lumbar vertebrae and the lower thoracic vertebrae of normal human bodies and it has a high degree of automation. This study is of great significance to the morphological study of the spine and the simulation of the scoliosis operation.

computed tomography image sequence; vertebra segmentation; area change regulation; three-dimensional reconstruction

TP 391

10.11996/JG.j.2095-302X.2019010040

A

2095-302X(2019)01-0040-06

2018-06-19；

2018-06-27

國家自然科學基金面上項目(61472466，61572167，61271123)；國家自然科學基金青年項目(61502136)

王 浩(1994-)，男，山東日照人，碩士研究生。主要研究方向為圖形圖像處理、軟件工程。E-mail：2013214020@mail.hfut.edu.cn

霍 星(1979-)，女，安徽合肥人，副教授，博士，碩士生導師。主要研究方向為圖形圖像處理等。E-mail：huoxing@hfut.edu.cn