基于空間變換與通用三維軟件的逆向參數(shù)化股骨建模技術(shù)研究

金帥普，王 玲，余 娟

（1.西安交通大學機械工程學院，陜西 西安 710049；2.西安郵電大學電子工程學院，陜西 西安 710121）

1 引言

目前，計算機數(shù)字化技術(shù)已廣泛應用于逆向工程中，特別是對于骨科臨床，已形成了新型數(shù)字骨科[1]。CT、X光和MRI等非接觸成像技術(shù)已成為獲得患者骨骼模型的主流方法。由于骨骼壞死、骨折等疾病嚴重影響病人的行動能力，常需要設計合適的骨科植入物、假肢等，在設計人工假體時，需要結(jié)合假體和病變骨骼模型進行動力學仿真或有限元仿真等一系列設計工作，需要一個易于編輯修改的三維骨骼模型。目前常通過CT掃描、Mimics圖像處理、Geomagic Studio逆向建模、Proe∕Ug三維建模[2]等方式進行三維骨骼的重建工作，步驟繁多，費時耗力，無法滿足這種要求。利用數(shù)學工具和三維設計軟件構(gòu)建一個可變形的參數(shù)化骨骼模型是一個具有挑戰(zhàn)性的工作。

現(xiàn)代醫(yī)學中，股骨頭壞死等病變是導致機體功能缺失的主要原因，為恢復生理功能，經(jīng)常采用假體置換的方法來構(gòu)造人工關(guān)節(jié)[3]。人工假體需要根據(jù)病人骨骼尺寸參數(shù)來定制，由于人的個性化差異，不同病人的骨骼具有不同特征參數(shù)。為此，提出以股骨為主要研究對象，構(gòu)建參數(shù)化股骨三維模型，使股骨關(guān)鍵參數(shù)可編輯的方法。該方法也可擴展到其他骨骼或具有復雜表面實體的參數(shù)化建模中。

文獻[4]利用曲面特征將股骨的網(wǎng)格模型分割為有醫(yī)學意義的曲面特征單元，設計特征參數(shù)及參數(shù)間約束關(guān)系，生成股骨參數(shù)化CAD曲面模型。文獻[5]使用Mimics、Imageware、Unigraphics等軟件對股骨點云數(shù)據(jù)進行處理，得到了較為精確的股骨三維模型。文獻[6]利用3D掃描的方法創(chuàng)建骨骼網(wǎng)格模型，該方法比2D掃描更有效。以上方法或無法精確重建骨骼模型的細節(jié)或步驟繁瑣，無法實現(xiàn)自動參數(shù)化。文獻[7]通過對點云數(shù)據(jù)進行簡單B樣條插值構(gòu)建整體曲面模型，未能將骨骼曲面參數(shù)表示。文獻[8]等依據(jù)95例國內(nèi)臨床計算機CT數(shù)據(jù)，通過網(wǎng)絡投影變換和統(tǒng)計學分析建立了能夠反映個體特征差異的股骨和脛骨的參數(shù)化模型，該方法無法實現(xiàn)對骨骼關(guān)鍵參數(shù)編輯。

綜合上述考慮，提出了一種基于Matlab和Creo的逆向參數(shù)化骨骼建模方法，如圖1所示。通過該方法，可獲得關(guān)鍵參數(shù)可編輯的骨骼三維模型，只需在Creo三維軟件中進行調(diào)整，無需多個軟件協(xié)同操作，方便快捷。通過實驗驗證，確定了該方法所建立骨骼三維模型的有效性。

圖1 逆向參數(shù)化骨骼建模流程圖Fig.1 Reversed Parametric Bone Modeling Flow Chart

2 參數(shù)化骨骼建模流程

股骨逆向參數(shù)化建模主要分為五步：（1）獲取點云數(shù)據(jù)；（2）點云數(shù)據(jù)的數(shù)學處理——空間變換[9]；（3）在通用三維軟件中進行實體建模；（4）添加控制參數(shù)；（5）驗證模型有效性?，F(xiàn)將其描述如下。

2.1 獲取點云數(shù)據(jù)

為獲取標準股骨點云數(shù)據(jù)，須先將標準股骨模型進行CT掃描，將掃描圖像放在Mimics軟件中進行圖像處理，根據(jù)灰度值不同[10]獲得圖像的點云數(shù)據(jù)。將點云數(shù)據(jù)導入到Creo軟件中，根據(jù)股骨表面曲率變化，分區(qū)域截取點云數(shù)據(jù)，曲率變化較小的區(qū)域截取較少截面的點云數(shù)據(jù)，曲率變化較大的區(qū)域截取較多截面的點云數(shù)據(jù)。

2.2 點云數(shù)據(jù)坐標變換

股骨點云數(shù)據(jù)，如圖2所示。股骨干建模坐標系易用Csys-1，股骨頸與股骨頭建模坐標系易用Csys-2，以這兩個典型坐標系之間的轉(zhuǎn)換關(guān)系為例進行說明，其他特征的局部坐標系不再贅述（例如大轉(zhuǎn)子區(qū)域局部坐標系）。為方便點云數(shù)據(jù)的獲取，通常以Csys-1作為基準坐標導出點云三維數(shù)據(jù)。當對股骨頸與股骨頭建模時，須將Csys-1下獲得的數(shù)據(jù)轉(zhuǎn)換到Cs ys-2下，以便后續(xù)三維軟件的建?？刂坪蛥?shù)化控制。本小節(jié)數(shù)學處理均在專業(yè)數(shù)學軟件Matlab中進行。

2.2.1 坐標變換原理

如圖2所示，將Cs ys-2中點M'(x2，y2，z2)變換到Cs ys-1中點M(x1，y1，z1)，其變換方程為：

式中：Csy s-1O Csys-2—Csys-2坐標原點在Cs ys-1中的坐標；—Csys-2變換到Csys-1的旋轉(zhuǎn)矩陣；由繞坐標軸的旋轉(zhuǎn)矩陣Rot(x，α)，Rot(y，β)，Rot(z，γ)按照一定順序相乘得到，其表示如下：

式中：c?=cos(?)，s?=sin(?)，α，β，γ—由Csys-2變換到Cs ys-1時繞x，y，z軸的旋轉(zhuǎn)角度。坐標變換示意圖，如圖3所示。

圖3 坐標變換示意Fig.3 Coordinate Transformation

2.2.2 坐標變換參數(shù)的確定

仍以Csys-1變換到Csys-2為例，參數(shù)Csys-1O Csys-2和α，β，γ確定如下。

（1）將股骨頸與股骨頭處等切成N份切片，將每個切片上的點云輪廓看作同一平面，運用最小二乘法[11]得到擬合平面。平面方程表達式為：

最小二乘法求解表達式中的a0，a1，a2，即：

（2）將切片點云數(shù)據(jù)映射到（1）中所求平面上。設投影點為，則向量平行于（1）中平面的法向量，則可得：

（3）對每3個不共線的投影點求圓心，取圓心平均值作為為該切片中心。設三點不共線，則線段的中垂面方程以及投影面方程聯(lián)立［12］：

解三個平面方程可得切片中心坐標C(x c，y c，z c)。

（4）由N個 切 片 可 得 到N個 中 心C k(x ck，y ck，z ck)，k=(1，2，3，…，N-1，N)，N個中心近似在一條空間直線上，形成股骨頸與股骨頭軸線。由最小二乘原理求取頸線L，設L為：

則需求取參數(shù)：x0，y0，m，n，化簡方程為：

則：

由最小二乘原理得：

解方程可得頸線L參數(shù)x0，y0，m，n。

（5）將Cs ys-1中z軸旋轉(zhuǎn)到與頸線L重合，先將其繞x軸旋轉(zhuǎn)α角度，再將其繞y軸旋轉(zhuǎn)β角度。α是頸線L在Csys-1中o yz面的投影L1與oz軸的夾角，β是頸線L與Csys-1中面o yz的夾角。頸線L的方向向量為的方向向量為，oz軸方向向量為，oyz的法向量為，則：

由于不繞z軸旋轉(zhuǎn)，所以γ=0。

（6）計算股骨頭球心。選取股骨頭3個不共面的點云數(shù)據(jù)B1(x1，y1，z1)，B2(x2，y2，z2)，B3(x3，y3，z3)，設股骨頭球心坐標為C(x c,y c,z c)，則，且C在頸軸線L上，可得：

解之可得股骨頭球心坐標C(x c，y c，z c)。取P組3個不共面點云數(shù)據(jù)，聯(lián)立P組方程組，解出P個球心坐標，取P組坐標平均值作為最終球心坐標C1(x c1，y c1，z c1)。取股骨頭球心點C1(x c1，y c1，z c1)作為Csys-2坐標原點，則Csy s-1O Csys-2=C1(x c1，y c1，z c1)。

由上述六步可最終確定坐標變換所需要的所有參數(shù)C sys-1O Cs ys-2和α，β，γ，參數(shù)確定示意圖，如圖4所示。將（2）中映射點云M'i(x'i，y'i，z'i)進行坐標轉(zhuǎn)換，部分轉(zhuǎn)換數(shù)據(jù)，如表1所示。

表1 坐標轉(zhuǎn)換數(shù)據(jù)Tab.1 Coordinate Transformation Data

圖4 參數(shù)確定Fig.4 Parameter Determination

2.2.3 笛卡爾坐標轉(zhuǎn)極坐標

考慮到通用三維軟件Creo中用角度來控制截面輪廓的生成，須將點云數(shù)據(jù)變換為(0～360)°范圍內(nèi)的極坐標數(shù)據(jù)。將前述坐標轉(zhuǎn)換后數(shù)據(jù)分別在各自坐標系統(tǒng)中進行極坐標轉(zhuǎn)換，其變換公式為：

式中：θ—點的極角；x、y—點的橫軸坐標值和縱軸坐標值；?—待確定的角度參數(shù)。由于θ的取值范圍是(0～360)°，而反正切三角函數(shù)得出的取值范圍為(-90～90)°，所以加上?使其值域為(0～360)°。?由點所在的象限來確定。當x≥0且y≥0時，?=0；當x<0且y>0時，?=90°；當x<0且y<0時，?=180°；當x>0且y<0時，?=270°。

2.3 建立三維模型

將點云坐標轉(zhuǎn)換到合適的坐標系統(tǒng)后，即可在三維軟件中進行建模。采用Creo三維參數(shù)化軟件進行股骨建模。

2.3.1 建立“圖形”特征與截面輪廓特征

在三維軟件Creo中，“圖形”特征可以看作用來控制“掃描”特征的函數(shù)?！皥D形”特征由樣條曲線產(chǎn)生，樣條曲線由極坐標點云數(shù)據(jù)控制。

（1）創(chuàng)建“圖形”特征。進入Creo軟件的“圖形”特征模塊，繪制任意樣條曲線，將樣條曲線首尾端點之間的距離設置為360，對應一個封閉輪廓的360°取值范圍。將該“圖形”特征命名為“dy-k”，如圖5中間圖形。

圖5 截面輪廓特征Fig.5 Profile Featur

（2）導入前述極坐標值。雙擊樣條曲線，導入前述第k切片極坐標值，用于控制樣條曲線的形狀。約束樣條曲線首尾端點處與豎直軸為“垂直”關(guān)系，使封閉曲線的0°點和360°點光滑連接。

（3）創(chuàng)建截面輪廓特征。以股骨干輪廓建模為例，首先建立Csys-1基準坐標系，該坐標系不需要參考，可將Creo軟中默認坐標系作為Csys-1。建立平行于xo y面偏移距離為z值大小基準平面。在此基準平面上繪制半徑為2.5mm，圓心為Csys-1原點的圓。進入“掃描”功能模塊，創(chuàng)建圖5右側(cè)圖形中“起點”處黑色線段，假設該線段尺寸變量名稱為sd4，則創(chuàng)建關(guān)系式：

其中，evalgraph函數(shù)用于關(guān)聯(lián)“圖形”特征與掃描截面，tra‐jpar是取值為[0，1]的變量。每一個θ對應一個ρ，旋轉(zhuǎn)360°即可掃描出第k個切片的截面輪廓。

2.3.2 建立“混合”特征

創(chuàng)建出每個切片處的輪廓后即可“混合”出三維實體。打開Creo軟件中的“混合”功能，選取前述生成的各個截面輪廓，將“混合曲面”類型選取為“平滑”，設置起始截面邊界條件為“相切”以便不同“混合”特征之間的平滑連接?！盎旌稀苯Ｊ疽鈭D，如圖6所示。

圖6 “混合”建模Fig.6 Solid Modeling

2.4 添加控制參數(shù)

在進行股骨動力學建模和有限元分析時，主要關(guān)注五個參數(shù)，分別是：股骨頭半徑R、頸干角A、偏心距E1、股骨干長度L、橫向縮放因子K，如圖7左側(cè)圖形。創(chuàng)建這些參數(shù)的過程如下：

（1）添加參數(shù)。打開Creo中“參數(shù)”選項，添加五個參數(shù)，根據(jù)標準股骨模型尺寸賦予初始參數(shù)值，添加參數(shù)說明。參數(shù)設置如圖7右側(cè)圖形。

（2）添加參數(shù)與尺寸的關(guān)系。Creo中每個尺寸都有一個變量名稱，將變量名稱與（1）中的參數(shù)建立關(guān)系即可用參數(shù)驅(qū)動模型。

①股骨干長度參數(shù)化。股骨干長度和每一層切片對應的坐標z有關(guān)，假設股骨頸軸線與股骨干軸線交點處切片高度為z0，則建立關(guān)系：

式中：sdi—第i層切片對應z坐標的尺寸變量名稱，即基準平面偏移距離。則改變L值整體股骨干程度均增大一定比例，最終長度變?yōu)長。

②頸干角參數(shù)化。結(jié)合3.2.2中（5）以及圖7左側(cè)圖形可知：A=180°-α，設初始模型中α=α0，則建立關(guān)系：

圖7 控制參數(shù)Fig.7 Control Parameter

式中：sd2—頸干角A的尺寸變量名稱。改變A值相當于在初始頸干角的基礎上增加改變量，最終使頸干較為改變后的A值。

③偏心距參數(shù)化。由圖2可知Csys-2相對于Cs ys-1在y方向的平移距離y c1即為偏心距，建立關(guān)系：

式中：sd1—偏心距尺寸變量名稱。改變E1值可改變股骨頭中心相對股骨干軸線偏移距離。

④股骨頭半徑參數(shù)化。建立關(guān)系：

式中：s d3—股骨頭半徑R的尺寸變量名稱。改變R值，則股骨頭半徑隨之改變。

⑤橫向縮放因子。建立關(guān)系：

式中：sd4—3.3.1中（3）掃描截面線的尺寸變量名。改變k值，則股骨橫向等比例縮放，從而控制股骨的橫向尺寸。

2.5 驗證模型有效性

通過改變參數(shù)進行實例測試和誤差分析，闡述了該方法及所構(gòu)建模型的準確性。

2.5.1 構(gòu)建實例

將參數(shù)化股骨模型按照表2中參數(shù)進行變換對比。參數(shù)變化前后三維模型對比圖，如圖8～圖10所示。

圖8 股骨模型參數(shù)改變前后對比圖Fig.8 Contrast Graphics Before and After Changes in Femoral Model Parameters

利用參數(shù)與尺寸之間的變化關(guān)系，實現(xiàn)股骨的關(guān)鍵參數(shù)變形。由圖8～圖10可以看出，股骨干長度的改變導致了大小轉(zhuǎn)子高度的改變，股骨頸軸線也隨參數(shù)A的變化而變化，股骨頭半徑在改變后明顯增大，偏心距同樣增大。

2.5.2 誤差分析

為驗證本模型的有效性，如圖9所示。使用所構(gòu)建的特征模型與一個按照傳統(tǒng)方式由CT數(shù)據(jù)構(gòu)建的股骨樣本參數(shù)，如表2中更改后參數(shù)所示。模型進行對比，并采用表面積和體積描述模型與樣本的形狀相似度量。

圖9 股骨模型與傳統(tǒng)方式模型體積與表面積對比圖Fig.9 The Comparison Between Volume and Surface Area

表2 參數(shù)對照表Tab.2 Parameter Comparison Table

所構(gòu)建的股骨參數(shù)化三維模型在絕大部分情況下可準確表示股骨解剖特征，如表3、圖9所示。比較所建模型與傳統(tǒng)方式所建模型可知[13]，表面積相對誤差為0.27%，體積相對誤差為0.03%，根據(jù)經(jīng)驗可知，這種誤差不會對假體設計中的動力學建模和有限元仿真造成較大影響。

表3 誤差統(tǒng)計表Tab.3 Error Statistics

中紅色部分為所建模型與傳統(tǒng)方法所建模型對比中Haus‐dorff距離大于1mm的區(qū)域[14]，如圖10所示。

從圖10中可以看出股骨頭部分誤差較大，這是因為參數(shù)化模型中頭部為半球模型，而股骨存在一個假體設計中忽略的股骨凹，其對動力學建模和有限元仿真影響不大。大轉(zhuǎn)子處誤差也較大，因為大小轉(zhuǎn)子存在多種結(jié)構(gòu)[15]，但其對小轉(zhuǎn)子部分骨折的接骨板設計影響不大。股骨干側(cè)面誤差較大是因為選取的切片密度太小，沒有很好地擬合股骨干的局部特征，增大切片密度可有效解決該問題。

圖10 Hausdorff距離大于1mm的區(qū)域Fig.10 Hausdorff Distance Greater than 1mm

3 結(jié)論

針對人體骨骼的復雜曲面模型，提出了一種基于特征參數(shù)化技術(shù)自動構(gòu)建股骨三維結(jié)構(gòu)的方法。與傳統(tǒng)獲得人體股骨模型方法相比，具有以下特點：（1）由最小二乘原理確定股骨切片點云數(shù)據(jù)所在平面；（2）由幾何關(guān)系和線性擬合確定股骨頭與股骨頸的軸線；（3）通過參數(shù)驅(qū)動尺寸實現(xiàn)股骨三維模型的自動化創(chuàng)建。該方法對其他復雜曲面實體模型的逆向建模工作有借鑒意義。

股骨三維模型的自動化創(chuàng)建仍有很多尚未解決的問題，目前只對股骨特征提取和關(guān)鍵參數(shù)的自動化進行了初步研究。下一步工作是對患者股骨細節(jié)特征進行研究，進一步提高所建三維模型的精度，便于指導骨科假體和手術(shù)植入物的設計。