裴欣欣，羅嗣頻，萬業(yè)達(dá)，王學(xué)民,3，周 鵬,3*(.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 30007；.天津醫(yī)院放射科，天津 300； 3.天津市生物醫(yī)學(xué)檢測技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室，天津 30007)

基于連續(xù)成像的頸椎運(yùn)動特征

裴欣欣1，羅嗣頻2，萬業(yè)達(dá)2，王學(xué)民1,3，周 鵬1,3*
(1.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072；2.天津醫(yī)院放射科，天津 300211； 3.天津市生物醫(yī)學(xué)檢測技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室，天津 300072)

目的 在獲取X線頸椎屈伸運(yùn)動連續(xù)圖像的基礎(chǔ)上探討頸椎的運(yùn)動特征。方法 對頸椎病患者及正常人頸椎過屈位至過伸位自發(fā)連續(xù)運(yùn)動的多幀矢狀位X線圖像進(jìn)行圖像預(yù)處理、配準(zhǔn)后，記錄單椎體的運(yùn)動軌跡，提取各椎體的特征點(diǎn)，并計算出活躍椎體的相對幾何學(xué)參數(shù)、運(yùn)動速率。結(jié)果 患者組與正常組C4、5、6椎體動態(tài)軌跡不同。與正常組比較，過屈位時患者組C4、5椎體左側(cè)邊緣曲線與基準(zhǔn)線的夾角較小(P均<0.05)?；颊呓MC4椎體的運(yùn)動存在不穩(wěn)定性，角度變化的波動性、速率改變的幅度及頻率均較大。且患者組相鄰單椎體相對位置的改變更小。結(jié)論 通過本研究圖像預(yù)處理、配準(zhǔn)、參數(shù)提取及結(jié)果分析方法，可真實(shí)反映頸椎病變帶來的功能學(xué)特征改變，為基于X線圖像的頸椎動態(tài)特征分析提供新思路。

頸椎；動態(tài)特征；運(yùn)動速率；單椎體運(yùn)動

傳統(tǒng)的頸椎病檢查多采用X線進(jìn)行形態(tài)學(xué)特征分析[1-2]，而頸椎病初期是功能性的，靜態(tài)觀察不能很好地反映其特征。以往通過有限元法[3]可追蹤到靜態(tài)圖像椎體邊緣，基于雙目立體視覺的頸椎運(yùn)動測量[4]、光學(xué)追蹤[5]和電磁追蹤[6-7]等均依賴于外部標(biāo)記運(yùn)動追蹤，所得數(shù)據(jù)誤差較大。Anderst等[8]將直徑1.6 mm的鈦珠植入研究對象的椎體和棘突內(nèi)以獲取椎體運(yùn)動特征，但方法具有侵入性。葉斌[9]通過攝取頸椎張口位、左右旋轉(zhuǎn)后張口位研究頸椎的動態(tài)運(yùn)動。但此類研究多采用不同位置的靜態(tài)矢狀位圖像，分別對圖像進(jìn)行測量分析，其動態(tài)特征不明顯。本研究通過訓(xùn)練受試者自主完成頸椎屈伸運(yùn)動，選取頸椎過屈位至過伸位自發(fā)連續(xù)運(yùn)動的多幀矢狀位X線圖像，經(jīng)預(yù)處理、配準(zhǔn)后，提取頸椎椎體特征點(diǎn)，通過與相鄰圖像的連續(xù)處理，描繪每一椎體的整體運(yùn)動的軌跡，計算椎體運(yùn)動速率，分析感興趣椎體動態(tài)功能特征和相鄰兩椎體在運(yùn)動過程中相對位置的改變，以期更真實(shí)地反映頸椎病變引起的頸椎動態(tài)功能的改變，進(jìn)一步加深對生理載荷下頸椎運(yùn)動狀況的認(rèn)識。

1 試驗(yàn)流程

在獲取頸椎屈伸運(yùn)動連續(xù)圖像的基礎(chǔ)上，對圖像進(jìn)行預(yù)處理，以增強(qiáng)頸椎的結(jié)構(gòu)，提取出每一椎骨的特征點(diǎn)。由于受試者進(jìn)行頸椎屈伸運(yùn)動時，C7椎體常被遮擋，X線圖像中難以區(qū)分，為清楚地了解各椎骨的相互關(guān)系，以C6椎體圖像進(jìn)行配準(zhǔn)處理。完成上述處理后，對相鄰的圖像進(jìn)行參數(shù)提取和計算，分析正常人群與頸椎病患者頸椎運(yùn)動速率的改變，并記錄單椎體的運(yùn)動軌跡和相對位置的變化，試驗(yàn)流程見圖1。

圖1 試驗(yàn)流程圖

2 圖像獲取

于2015年10月—12月征集21名受試者納入研究，包括臨床診斷為頸椎病患者14例(患者組)及健康志愿者7名(正常組)?；颊呓M中，男4例，女10例，年齡22～54歲，中位年齡37.0歲；正常組中，男2例，女5例，年齡23～45歲，中位年齡34.0歲。所有受試者對試驗(yàn)風(fēng)險均知情同意，自愿參與試驗(yàn)。

采用島津Sonial vision SafireⅡ平板數(shù)字化X線多功能透視攝影系統(tǒng)進(jìn)行頸椎動態(tài)檢測。設(shè)置管電流為2 mA，管電壓為60 kV，X球管與受試者頸部皮膚間的距離為1.2 m，調(diào)整成像中心于受試者頸椎C4椎體。為觀察頸椎的運(yùn)動狀態(tài)，X線攝片時囑受試者側(cè)位站立，雙手叉腰，肘向后靠攏，首先進(jìn)行低頭至過屈位，而后在15 s內(nèi)從過屈位勻速運(yùn)動至過伸位。將X線機(jī)置于“準(zhǔn)備”態(tài)，在受試者頸椎以勻速從過屈位運(yùn)動至過伸位的過程中進(jìn)行X線曝光，圖像采集速率為6幀/秒，共采集15 s[10]。所獲圖像為1179×782像素，將圖像按順序編號命名，取每一受試者8幀不同位置圖像，見圖2。

3 圖像預(yù)處理及配準(zhǔn)

3.1選取ROI 為方便后續(xù)的參數(shù)提取，對所獲取圖像進(jìn)行預(yù)處理，在Matlab中讀取已保存的圖像，選取涵蓋頸椎過屈位至過伸位的全部位置區(qū)域的最小圖像作為ROI，記為W，大小為u×v像素，其中u為橫向像素個數(shù)，v為縱向像素個數(shù)，對每一受試者選取的ROI圖像大小不同，但需使每個位置的頸椎均在ROI內(nèi)，圖3所示即為從圖2的圖像中選取的ROI。

3.2圖像增強(qiáng) 為提高圖像對比度，調(diào)整圖像的灰度直方圖分布，使頸椎部分更為突出且減少其周圍組織和噪聲的影響，根據(jù)圖像自身的特點(diǎn)，選取基于形態(tài)學(xué)的頂帽、底帽變換圖像增強(qiáng)算法和基于Retinex色彩理論[11]的圖像增強(qiáng)算法進(jìn)行雙增強(qiáng)。

頂帽、底帽變換的一個重要用途是校正不均勻光照[12]，在多張圖像的處理中，該算法對背景光照進(jìn)行了校正,見圖3。Retinex色彩理論是基于人類視覺基礎(chǔ)的理論，該算法具有銳化圖像、顏色恒常性、動態(tài)壓縮范圍大、可對不同圖像進(jìn)行自適應(yīng)增強(qiáng)等優(yōu)點(diǎn)[11]，但對噪聲較敏感[12]。聯(lián)合兩種方法可使圖像增強(qiáng)效果更好(圖4)，以便頸椎邊緣的完整提取。

3.3邊緣提取 為提取頸椎各椎體左側(cè)邊緣曲線得到頸椎運(yùn)動軌跡，從而觀察整體運(yùn)動趨勢，采用性能優(yōu)良Canny算子[13]初步提取邊緣，見圖5。

圖2 受試者頸椎過屈位(A～D)運(yùn)動至過伸位(E～H)X線圖像

圖3 頂帽、底帽變換圖像增強(qiáng)算法后圖像 圖4 Retinex色彩理論增強(qiáng)后圖像 圖5 采用Canny算子初步提取頸椎邊緣曲線 圖6 修正后提取的椎體左側(cè)邊緣曲線

在第1幀圖像的頸椎椎體邊緣選取m個種子點(diǎn)，并將各點(diǎn)坐標(biāo)保存于m×2的矩陣A(x,y)內(nèi)。在第2幀圖像內(nèi)遍歷像素點(diǎn)，在m個種子點(diǎn)的增長區(qū)域內(nèi)尋找像素值為1的點(diǎn)，并將此點(diǎn)的坐標(biāo)保存在與圖像大小相等的矩陣內(nèi)。種子點(diǎn)(xij,yij)的增長區(qū)域?yàn)椋?/p>

橫坐標(biāo)：(xij-a1,xij+a2),i=1,2,3,…,n;j=1,2,3,…,m

縱坐標(biāo)：(yij-b1,yij+b2),i=1,2,3,…,n;j=1,2,3,…,m

其中，a1,a2,b1,b2由多次實(shí)驗(yàn)得出，范圍為3～7，不同受試者的參數(shù)選擇不同。n為圖像的幀數(shù)，以包括過屈位和過伸位為邊界，一般為80～90幀。初次完成后，對提取的邊緣作出修正，此時可較光滑地提取頸椎椎體邊緣的左側(cè)曲線，見圖6。左側(cè)曲線的像素點(diǎn)即為新的種子點(diǎn)，將其保存于矩陣A內(nèi)，作為下一幀圖像的遍歷的種子點(diǎn)，各幀圖像均以前一幀圖像頸椎左側(cè)邊緣曲線為種子點(diǎn)進(jìn)行遍歷，各幀圖像依次進(jìn)行，并將各幀所得矩陣A內(nèi)的種子點(diǎn)坐標(biāo)均保存于數(shù)組C內(nèi)，數(shù)組C大小與W相同。

3.4邊緣曲線配準(zhǔn) 為觀察頸椎整體運(yùn)動的范圍和趨勢變化，減少在拍攝時受試者肩部位置改變的影響，選取第1幀圖像頸椎的C6椎體左下側(cè)頂點(diǎn)為基點(diǎn)進(jìn)行配準(zhǔn)，在提取完頸椎椎體左側(cè)曲線后，此點(diǎn)即為縱坐標(biāo)y值的最大點(diǎn)，坐標(biāo)記為M1(X1,Y1)，第2幀圖像相同位置的點(diǎn)坐標(biāo)記為M2(X2,Y2)，以此類推，最后一幀相同位置點(diǎn)坐標(biāo)記為Mn(Xn,Yn)。

配準(zhǔn)距離差為：Di(X,Y)=Mi(Xi,Yi)-M1(X1,Y1)；i=2,3,4…，n

配準(zhǔn)過程為：Bi(xij,yij)=Ai(xij,yij)-Di(X,Y)；i=2,3,4,…,n;j=1,2,3,…,m

表1 患者組與正常組C4、5椎體過屈位、過伸位的角度比較(°)

圖7 C4、5椎體左邊緣曲線 A.正常人； B.頸椎病患者

此時，頸椎椎體的左側(cè)曲線全部以C6椎體的左下側(cè)頂點(diǎn)為基準(zhǔn)點(diǎn)配準(zhǔn)，觀察頸椎整體的運(yùn)動趨勢、頸椎的生理曲線狀況及頸椎后伸位的最大位置所在，見圖7。

4 參數(shù)提取與結(jié)果分析

將各受試者的圖像歸一化為86幀，以C4、5椎體為感興趣節(jié)段，分別以其第1幀左側(cè)邊緣曲線為起點(diǎn)選擇種子點(diǎn)，用上述方法處理(圖7)，紅色代表配準(zhǔn)前所在位置，其他不同顏色代表各幀圖像各椎體左側(cè)邊緣所在位置。

患者組與正常組C4、5、6椎體動態(tài)軌跡不同，患者組頸椎生理彎曲消失，且C4、5椎體間隙較小，且在頸椎屈伸運(yùn)動過程中椎體間隙大小發(fā)生明顯變化。

記錄單椎體運(yùn)動，選取數(shù)組C內(nèi)第5～10行所有數(shù)據(jù)，按列求取平均值的相對固定點(diǎn)坐標(biāo)Ni(xi,yi),

i=2,3,…,n。為觀察單椎體的運(yùn)動趨勢，以配準(zhǔn)點(diǎn)的橫坐標(biāo)作基準(zhǔn)線，計算椎體左側(cè)邊緣曲線與基準(zhǔn)線的夾角：

i=2,3,4,…,n

患者組及正常組的C4與C5椎體的夾角及其在運(yùn)動過程中的變化見圖8、9。

以C4椎體分析頸椎的運(yùn)動速率，其角度變化見圖8B；對離散數(shù)據(jù)進(jìn)行差分運(yùn)算，其速率變化見圖10。由于頸椎病患者不能很好地控制頸椎的運(yùn)動從而產(chǎn)生被動的速率變化，表現(xiàn)為運(yùn)動的不穩(wěn)定性?；颊呓M角度變化的波動性大，其速率改變的幅度及改變頻率均較大(圖10)。患者組與正常組頸椎運(yùn)動速率方差對比見圖11，患者組方差為0.462，正常組為0.331。

觀察單椎體運(yùn)動軌跡(圖8、9)，分析患者組與正常組C4、5椎體左側(cè)邊緣曲線與基準(zhǔn)線的夾角，對患者組和正常組C4、5椎體過屈位、過伸位的角度采用SPSS軟件進(jìn)行兩獨(dú)立樣本t檢驗(yàn)結(jié)果見表1?；颊呓M與正常組C4、5椎體在過屈位時角度有差異有統(tǒng)計學(xué)意義(P均<0.05)，而在過伸位時差異無統(tǒng)計學(xué)意義(P均>0.05)。

分析相鄰單椎體相對位置的改變時，將C4、5椎體角度的變化相減，計算C4、5椎體相對位置變化的數(shù)據(jù)并擬合出模板?；颊呓M及正常組C4、5椎體相對位置變化數(shù)據(jù)比較見圖12，紅色為正常組模板數(shù)據(jù)，藍(lán)色為患者組的數(shù)據(jù)，可見患者組椎體相對位置的改變較小，曲線較正常組平緩且最低點(diǎn)絕對值較小。

圖8 C4椎體運(yùn)動軌跡 A.C4椎體夾角示意圖； B.C4椎體角度變化對比圖 圖9 C5椎體運(yùn)動軌跡 A.C5椎體夾角示意圖； B.C5椎體角度變化對比圖

圖10 正常人與患者椎體角度變化速率對比圖 圖11 正常人與患者椎體角度變化速率方差對比圖 圖12 正常人與患者C4與C5椎體相對位置變化

5 結(jié)語

本研究在對圖像進(jìn)行預(yù)處理的基礎(chǔ)上對頸椎椎體左側(cè)邊緣進(jìn)行提取分析，以配準(zhǔn)點(diǎn)橫坐標(biāo)為基準(zhǔn)線，觀察C4、5椎體在自發(fā)性頸椎屈伸運(yùn)動過程中椎體左側(cè)邊緣曲線與基準(zhǔn)線夾角的變化趨勢，分析頸椎運(yùn)動速率，發(fā)現(xiàn)頸椎病患者頸椎運(yùn)動過程中速率呈跳躍性變化。此外，本研究以C4、5椎體左側(cè)邊緣與基準(zhǔn)線的夾角的差值反映單椎體在運(yùn)動過程中相對位置的改變，基于連續(xù)圖像的動態(tài)特征分析，為單椎體運(yùn)動的研究提供更多數(shù)據(jù)，也為頸椎病診斷提供更多參考信息。但本研究的樣本量較少，且因存在個體差異，不同患者C4、5椎體在運(yùn)動過程中相對位置變化曲線不同,該差異還有待進(jìn)一步深入研究。

[1] Nagamoto Y, Ishii T, Sakaura H, et al. In vivo three-dimensional kinematics of the cervical spine during head rotation in patients with cervical spondylosis. Spine (Phila Pa 1976), 2011,36(10):778-783.

[2] Nagamoto Y, Iwasaki M, Sugiura T, et al. In vivo 3D kinematic changes in the cervical spine after laminoplasty for cervical spondylotic myelopathy. J Neurosurg Spine, 2014,21(3):417-424.

[3] Xu X, Hao HW, Yin XC, et al. Automatic segmentation of cervical vertebrae in X-ray images. International Joint Conference on Neural Networks, 2012,20:1-8.

[4] 鄭偉龍,王志勇,王世斌,等.頸椎三維運(yùn)動測量及生物力學(xué)特性研究//北方七省市區(qū)力學(xué)學(xué)術(shù)會議，2014：70-71.

[5] Lee JH, Kim JS, Lee JH, et al. Comparison of cervical kinematics between patients with cervical artificial disc replacement and anterior cervical discectomy and fusion for cervical disc herniation. Spine J, 2014,14(7):1199-1204.

[6] Bell KM, Bechara BP, Hartman RA, et al. Influence of number of operated levels and postoperative time on active range of motion following anterior cervical decompression and fusion procedures. Spine (Phila Pa 1976), 2011,36(4):263-268.

[7] Bechara BP, Bell KM, Hartman RA, et al. In vivo analysis of cervical range of motion after 4- and 5-level subaxial cervical spine fusion. Spine (Phila Pa 1976), 2012,37(1):E23-E29.

[8] Anderst WJ, Vaidya R, Tashman S. A technique to measure three-dimensional in vivo rotation of fused and adjacent lumbar vertebrae. Spine J, 2008,8(6):991-997.

[9] 葉斌.動態(tài)X線片在頸椎病早期診斷中的作用.現(xiàn)代診斷與治療,2013,24(5):1130-1131.

[10] Bifulco P, Cesarelli M, Romano M, et al. Measurement of intervertebral cervical motion by means of dynamic X-ray image processing and data interpolation. Int J Biomed Imaging, 2013,2013:152920.

[11] 陳琛,張建州.改進(jìn)的多尺度Retinex醫(yī)學(xué)X射線圖像增強(qiáng)算法.計算機(jī)工程與應(yīng)用,2015,51(9):191-195.

[12] 馬超玉.光照不均勻條件下圖像增強(qiáng)算法研究.長春:長春理工大學(xué), 2014:4-13.

[13] Wang XJ, Liu XM, Yong G. Image edge detection algorithm based on improved canny operator. Computer Engineering, 2012,38(14):168-172.

Motion characteristics of cervical spine based on continuous images

PEIXinxin1，LUOSipin2，WANYeda2，WANGXuemin1,3，ZHOUPeng1,3*
(1.SchoolofPrecisionInstrumentandOpto-ElectronicsEngineeringofTianjinUniversity,Tianjin300072,China; 2.DepartmentofRadiology,TianjinHospital,Tianjin300211,China; 3.TianjinKeyLaboratoryofBiomedicalTestingTechnologyandInstruments,Tianjin300072,China)

Objective To explore the motion features of cervical spine based on continuous X-ray images. Methods The cervical spontaneous continuous multi-frame sagittal images from flexion to extention positions were selected from cervical spondylosis patients (patients group) and healthy adult (normal group). After preprocessing and registration, the trajectory of single vertebral body were recorded, and the feature points of each vertebral body were extracted. Meanwhile， the relative geometry parameters and movement rate were calculated. Results The motion trajectory of patients' cervical spine C4—6 were different from healthy people. The angles between the left edge curve of the vertebral body (C4 and C5 vertebra) and the baseline of patients group were smaller than those of normal group in flexion position (allP<0.05). There were instability in the movement of C4 vertebral body in patient group, and the volatility of the angle change, the rate of change and the frequency were larger. The relative position change of the adjacent single vertebral body in the patient group are smaller. Conclusion Through the preprocessing, registration, parameter extraction and result analysis, the changes of functional features in cervical spondylosis patients are truly reflected. And it also provides a new idea for dynamic analysis of cervical vertebrae based on X-ray images.

Cervical vertebral; Dynamic characteristics; Movement speed; Single vertebral body motion

裴欣欣(1991—)，女，河北邢臺人，在讀碩士。研究方向：醫(yī)學(xué)圖像處理。E-mail: peixx@tju.edu.cn

周鵬，天津大學(xué)精密儀器與光電子工程學(xué)院，300072； 天津市生物醫(yī)學(xué)檢測技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室， 300072。E-mail: zpzp@tju.edu.cn

2016-10-17

2017-04-10

10.13929/j.1003-3289.201610065

TP391.5

A

1003-3289(2017)07-1090-05