三維橈骨成角楔形截骨術前自動規劃算法

石志良，廖詩旗，甘梓博，祝少博

（1.武漢理工大學 機電工程學院，武漢 430070；2.武漢大學中南醫院 骨科，武漢 430071）

0 引言

成角畸形是前臂骨畸形中一種常見的畸形類型。前臂骨折發生成角畸形時，改變了前臂固有的生理曲度，破壞了“旋前弓”和“旋后弓”，出現旋轉過程中骨性阻擋或因成角畸形使骨間膜緊張，因而限制了前臂旋轉活動［1］。目前，治療成人前臂骨折的金標準是解剖復位、穩定的鋼板固定和保護周圍的血液供應［2］。為了達到精確治療的目的，使用截骨矯形術可以有效矯正畸形、緩解疼痛、改善肢體功能和外觀。然而，傳統手動截骨矯形方法操作繁瑣，矯正準確率低，因此，以自動化方式生成術前規劃方案是截骨矯形手術的發展趨勢。

手工交互截骨矯形分為二維手工截骨和三維（Three-Dimensional，3D）手工截骨。傳統的二維（Two-Dimensional，2D）手術計劃一般是利用前后視圖、側位視圖的正交X 光片，由醫生進行人工軸線繪制和角度測量，以確定畸形部位、截骨方式和截骨位置［3］。文獻［4］中通過疊加兩側射線照片計算兩個平面中最大變形面和真實畸變角度作為術前計劃。這種依賴二維影像進行人工矯形的方法無法完全呈現三維解剖結構的真實情況，因此在畸形診斷上存在較大的限制。三維手工截骨矯形方法是在畸形骨與復位模板的近端或遠端配準的基礎上，對比畸形骨與健康骨的差異部分，根據醫生臨床經驗選擇合適位置進行截骨。文獻［5-8］中提出了使用三維計算機輔助手工截骨矯形方法，這些方法通過三維重建技術使截骨矯正過程可視化，使醫生能夠更直觀地了解規劃過程，從而提升治療效果和效率。但是，手工交互操作依賴于臨床參數和規劃者的經驗，很難精準定位截骨位置，反復實驗會浪費大量臨床成本。

針對現有手動矯形方法存在的缺陷，研究人員提出利用幾何計算與優化計算的方法自動生成術前規劃方案。幾何計算的方法通過畸形骨與對側健康骨的遠、近端配準變換矩陣解耦螺旋軸獲取旋轉角度，以解剖軸線的交點定義截骨位置進行截骨矯形。文獻［9］中提出了一種基于幾何計算的三維斜向單切旋轉截骨矯形方法，目的是確定最佳截骨位置和截骨面平移量恢復旋轉對齊；文獻［10-11］中深入探討了適用于不同類型長骨的三維術前自動截骨矯形規劃方法，這些方法能夠自動進行截骨位置和矯正角度的幾何計算，同時考慮了臨床上的限制因素。幾何方法通常依賴于給定的幾何規則和參數進行計算，對于復雜的非線性問題可能無法提供最優解，因此不能給出精確可行的術前規劃方案?；趦灮姆椒ǜ鶕R床約束和臨床目標，對優化算法中的參數及目標區域進行約束，建立基于臨床目標的數學模型，實現術前方案的自動生成。優化方法的優勢是適用于各種問題，計算結果可以根據給定的目標函數和約束條件進行調整和改進，以找到最優解或接近最優解的解決方案。針對橈骨旋轉截骨術，文獻［12］中將畸形骨遠端質心作為復位目標，使用單純下降法對截骨位置與復位角度優化計算。文獻［13］中采用斜向雙切旋轉截骨術，通過對骨長度、截面對齊等使用優化算法，將遠端、中部和近端骨段與目標骨自動對齊矯正復雜的骨骼變形。這類方法提供了一種精確、簡便且臨床可用的術前規劃方案，但均針對橈骨旋轉成角混合畸形，不適用于成角畸形。針對橈骨旋轉、成角楔形截骨術，文獻［14-15］中提出了一種橈骨自動截骨矯形的算法框架，矯正復位取得了較好的效果；但方法需要同時將截骨切平面、骨碎片復位、固定板的位置、螺釘的方向和位置作為優化目標，導致優化方法涉及的參數較多，計算成本急劇增加。

此外，在計算機輔助三維規劃手術中，健康對側常被用作規劃和評估的參考［16-17］。但是，大部分人使用慣用手更頻繁，造成了上肢不對稱。因此，文獻［11］中提出人體雙側前臂長度差異可以通過對模型縮放進行補償；文獻［18］中提出應用線性回歸模型描述橈骨和尺骨之間的雙邊差異，以糾正左右臂之間的長度差異；但這類用于補償前臂長度差異的方法尚未應用于橈骨截骨矯形的自動規劃問題。

在筆者的前期工作中，針對上肢畸形，提出了一種多階段、多目標的加權優化方法，以解決截骨面、復位對齊、固定板和螺釘空間位置等多個問題［19］；對于上肢成角旋轉混合畸形，利用遺傳算法對上肢的復位對齊和骨接觸面重疊率進行優化，實現了單刀斜切截骨矯形［20］；對于上肢成角畸形，采用擬牛頓法，通過確定截骨面和旋轉角實現楔形截骨矯形［21］；擬牛頓法具有較好的收斂性，但存在計算復雜、易受初始點影響和易陷入局部極小值等問題。因此，本文采用內點算法，提出了一種解決上肢成角畸形問題的優化算法。

基于上述問題，針對橈骨成角畸形，本文提出一種橈骨楔形截骨矯形術前自動規劃方法。首先，利用對側鏡像骨模型補償差異后的模型作為矯形復位模板，通過迭代最近點（Iterative Closest Point，ICP）算法對橈骨近端和遠端關節進行配準，自動計算畸形區域；其次，通過對橈骨遠端關節的解剖區域進行加權配準，確定三維畸形軸的方向；隨后，調整畸形骨使它與世界坐標系Y軸平行，并將畸形區域映射至XOZ平面，其中的畸形輪廓曲線采用三次樣條插值法確定，從而確定復位旋轉軸方位；最后，將關節解剖區域對齊作為截骨矯形的優化目標，以畸形區域為截骨平面約束、復位角的旋轉范圍為復位角約束，通過單目標優化算法進行迭代求解，自動生成最優的橈骨楔形截骨矯形術前規劃方案。

1 楔形截骨矯形術前自動規劃算法

1.1 算法流程

三維橈骨成角楔形截骨術前自動規劃算法的總流程如圖1 所示，具體步驟如下：

圖1 算法流程Fig.1 Algorithm flow

1）將CT 數據進行三維重建，生成畸形和對側健康橈、尺骨模型；

2）矢狀面作為對稱面，將對側健康的骨模型進行鏡像，隨后將鏡像得到的健康模型用作畸形骨矯形復位的模板；

3）采用迭代最近點配準算法對齊兩側橈骨和尺骨，計算兩側差異量并補償，重建新的矯形復位模板作為參考骨；

4）計算畸形橈骨的畸形區域范圍；

5）采用基于關節解剖區域權重的配準方法，計算橈骨三維畸形軸方向；

6）畸形骨變換至與世界坐標系Y軸平行，根據畸形區域范圍計算畸形輪廓在XOZ平面的曲線及復位旋轉軸方位；

7）畸形輪廓曲線上任取一點構建復位旋轉軸，將畸形區域范圍和復位旋轉角度范圍作為控制變量約束，以橈骨遠端解剖區域配準均方根誤差（Root Mean Square Error，RMSE）最小為目標設計目標函數，構建單目標優化模型；

8）優化迭代計算出最優的截骨平面和復位旋轉角度，自動生成橈骨楔形截骨矯形術前規劃方案。

1.2 補償雙側差異

1）通過CT 數據創建畸形和健康側橈骨和尺骨的三維三角形表面模型。

2）將三維世界坐標系中的XOZ平面作為對稱平面，獲取健康模型中每一個坐標點關于XOZ平面的對稱點，根據點云數據生成對側健康前臂的鏡像模型。

3）采用雙側補償的方法［18］糾正左右臂橈尺骨之間的長度差，提高術前規劃的復位精度。具體步驟如下：

a）計算橈骨和尺骨模型的質心，并確定三條相互垂直的主軸線通過這兩個質心。其中具有最小轉動慣量的骨軸即為重力軸，對應橈骨的主旋轉軸即橈骨長軸LR，對應尺骨的主旋轉軸即尺骨長軸LU，如圖2 所示。

圖2 補償雙側差異Fig.2 Compensation for bilateral differences

b）橈骨近端和遠端配準對齊得到矩陣Mpr和Mdr，由這兩個矩陣可獲得修正矩陣使橈骨遠端相對于近端重新定位。同時，尺骨近端和遠端配準對齊得到矩陣Mpu和Mdu，由這兩個矩陣可獲得補償矩陣

c）由補償矩陣Ma分解出平移矩陣，通過平移矩陣獲取變換過程中沿XYZ三個軸的位移變換量，計算出沿重力軸的位移量ΔZUlna。

其中：C=0.98 是比例常數；ΔZRadius是沿橈骨軸所需的補償半徑；ΔZUlna是沿骨軸的位移量。

1.3 畸形區域

在優化過程中計算截骨位置時，采用一種畸形區域的自動計算方法［13］能夠明確畸形骨偏離參考骨的具體位置，縮小搜索空間，同時減少計算量和時間成本。

1）利用ICP 算法對畸形及參考骨進行近、遠端配準，計算兩端的配準點集；

2）沿畸形骨模型的骨長軸將配準后的畸形骨模型和參考骨模型的點集以合適步長進行空間離散處理；

3）每一個步長內包含畸形骨模型和參考骨模型的點集，計算并搜尋畸形骨模型點集中每一點與參考骨模型點集中歐氏距離最小的點，通過計算每一個包圍盒中的兩個點集間的RMSE，衡量在此位置畸形骨與對側鏡像骨的偏離程度；

其中：PSD表示盒內畸形骨模型點集；PSN表示盒內參考骨模型點集；|PSD|表示集合PSD的基數。自動計算出的畸形函數圖像如圖3 所示。

圖3 骨畸形函數Fig.3 Bone deformity function

4）根據外科醫生建議，以偏離RMSE 最小值15%為畸形區域，計算畸形骨的畸形區域如圖4 中所示。

圖4 骨畸形區域Fig.4 Area of bone deformity

1.4 三維畸形軸

1.4.1 基于權重的關節局部區域配準

在臨床上，橈骨遠端由于與腕關節及尺骨遠端關節相連，它的關節面的對齊比長骨區域的對齊更重要。因此，采用基于權重的關節解剖區域配準方法［13］對畸形骨與參考骨遠端的標記區域進行精準復位對齊，以減小粗略對齊產生的誤差。

1）在遠端關節解剖區域選取7 個解剖點，并對7 個解剖點確定的解剖區域設置不同權重，權重取值范圍為［0，1］，具體分配情況見表1。

表1 解剖點序號對應解剖位置名稱及權重值Tab.1 Anatomic point serial number and corresponding anatomic position name and weight value

2）采用KDtree 算法［22］在7 個解剖點鄰域內搜索50 個點。ICP 配準的評價指標［23］：待配準點集與參考點集之間經過矩陣變換后，點集數據之間滿足某種度量準則下的最優匹配，即平均距離差值最小。因此，提出一種基于區域權重的點云配準誤差評價指標如下：其中：K表示橈骨遠端解剖區域的總數；N=50 為最鄰近點對個數；wi表示各個解剖區域的權重值；R表示從畸形骨解剖點集變換至參考骨解剖點集的旋轉矩陣；pij表示畸形骨解剖區域i中的第j個最近點，qij表示參考骨鏡像模型解剖區域i中的第j個最近點。

1.4.2 三維畸形軸方向

1）畸形骨和參考骨遠端通過關節整體區域配準，獲得變換矩陣Md；畸形骨和參考骨近端通過基于權重的關節局部區域配準，獲得變換矩陣Mp，由兩個變換矩陣產生校正矩陣

2）Mc矩陣包括平移和旋轉矩陣，從其中提取出不包含平移的矩陣R。

3）由矩陣R分解出3D 畸形軸的方向向量k［9］，如圖5 所示，計算方法如下：

圖5 3D畸形軸Fig.5 3D axis of malformation

1.5 畸形輪廓曲線和復位旋轉軸方位

1.5.1 旋轉軸方向向量

確定畸形軸的方向，構建同向的單位向量并進行變換，具體步驟如下：

1）以世界坐標系原點為起點，構建與三維畸形軸方向相同的單位向量V，并將單位向量V變換至與世界坐標系Y軸平行，即V=(0，1，0)；

2）計算方向向量V變換到與世界坐標系Y軸平行的變換矩陣MY；

3）對畸形骨及參考骨作MY變換，使其與Y軸平行，將畸形骨的畸形區域輪廓投影到XOZ平面，通過三次樣條插值求解畸形輪廓曲線函數表達式。

1.5.2 畸形輪廓曲線

通過投影到XOZ面的畸形區域輪廓范圍計算截骨平面的最近位置zmin和最遠位置zmax。對畸形區域分段三次樣條插值，其中兩個端點為z0=zmin，zn=zmax，選取合適步長將區間[zmin，zmax]分成n個區間，畸形輪廓線在X軸上的極值點xmin、xmax作為旋轉點的備選點，一側控制開放楔形截骨，一側控制閉合楔形截骨。計算過程如下：

在每一個小區間(zi-1，zi)(i=1，2，…，n)內，S(z)分別為三次多項式函數：

將數據節點和指定的首尾端點條件代入矩陣方程求解可得樣條曲線系數ai、bi、ci、di，即可計算出曲線方程S1(z)，曲線方程S2(z)?；屋喞€圖如圖6 所示。

圖6 畸形輪廓曲線Fig.6 Deformity contour curve

1.5.3 復位旋轉軸方位

以畸形輪廓曲線上任意一點(S(z)，0，z)作為復位旋轉軸位置，以旋轉軸方向V=(0，1，0)作為復位旋轉軸方向，確定復位旋轉軸方位，其中S(z) ∈[zmin，zmax]。

1.6 優化方法

采用內點算法將懲罰函數定義在可行域內，并從可行域內某一初始點出發，在可行域內進行迭代。該算法的優勢在于通過構造懲罰函數，將原有不等式約束優化問題轉化為一個在定義域內的無約束非線性規劃優化問題，結合牛頓法和基于信賴域法的共軛梯度方法進行求解。手動選取截骨位置及復位角計算繁瑣，且找到最優方案較為復雜，自動計算的截骨位置不佳無法確定截骨方式，因此通過內點算法進行最優解計算。以橈骨遠端解剖區域配準RMSE 最小為復位目標，畸形區域作為初始點的可行區域，輪廓曲線上任一點作為旋轉軸位置參數，橈骨遠端關節繞復位旋轉軸旋轉的角度為復位角參數，通過建立基于截骨平面位置和復位角變化的控制變量約束設計目標函數，利用邊界和約束條件控制算法的搜索空間，優化迭代計算出最優的截骨位置參數和復位角參數，達到最優的遠端關節解剖區域配準精度，自動生成橈骨楔形截骨矯形術前規劃。

1.6.1 變量約束

變量約束包括截骨位置約束和復位旋轉軸旋轉的角度約束。

截骨位置約束是指在截骨過程中截骨平面可變化的范圍。截骨平面決定截骨的大小和位置以及截骨的方向，不同的截骨位置會導致復位誤差不同。優化開始前，通過骨畸形自動診斷確定截骨平面大致范圍，即為截骨平面位置約束。

復位旋轉軸旋轉的角度約束是遠端骨段繞復位旋轉軸旋轉至與參考骨模型遠端的解剖區域對齊，即為復位旋轉軸旋轉的角度約束。旋轉的角度θ變化范圍為(-45°，45°)。

1.6.2 目標函數

橈骨成角畸形截骨術的術前優化目標為橈骨遠端的關節局部區域配準RMSE，用于測量配準精度，精細控制復位對齊。具體描述如下：

橈骨遠端關節解剖區域的配準精度：

其中：K表示解剖特征區域總數；wi表示權重，分配比例如表1；p表示畸形骨解剖區域的點，q表示重建模板解剖區域的點；f是畸形骨標志區域所有點pj與復位模板標志區域的點qj之間的歐氏距離的加權RMSE 值；R是將畸形骨的遠端骨段旋轉到正確解剖位置的變換矩陣。R的具體計算如下：

其中：T1表示將復位旋轉軸平移到坐標原點的平移矩陣，T2表示T1的逆變換，R(θ)表示將復位旋轉軸繞世界坐標系Y軸旋轉θ角度的旋轉矩陣。

1.6.3 內點算法

將實際問題轉化為算法求解的形式，具體步驟如下：

1）建立數學模型。式（13）為目標函數，式（14）為不等式約束集合：

2）引入松弛因子，構造對數障礙函數作為懲罰函數，生成等價于原模型的優化問題，并將不等式約束問題轉化為等式約束問題。

其中：r()k(k=0，1，2，…)為懲罰因子，是遞減的正數序列；r(k-1)=r(k)c，c為遞減系數；si≥0(i=1，2，3，4)為松弛因子。

3）在定義域內任取初始點x(0)，初始懲罰因子r(0)>0，0.1

4）構造拉格朗日函數，處理等式約束。

其中：Fr(z，θ，s)是目標函數的逼近問題；λj(j=1，2，3，4)為拉格朗日乘數；H為拉格朗日函數的Hessian 矩陣。

5）利用迭代法求解庫恩塔克條件（式（17）），獲得最優點。

求解過程中，當拉格朗日函數的Hessian 矩陣正定時，目標函數在迭代點附近鄰域內是嚴格凸函數，算法采用牛頓法。當Hessian 矩陣半正定時，目標函數在迭代點附近鄰域內是凹函數，牛頓法失效，算法轉為基于信賴域的共軛梯度法求解。逼近問題（15）在迭代附近不是局部凸時，直接采用基于信賴域的共軛梯度法。

6）若‖x?(r(k))-x?(r(k-1)) ‖≤ε或達到迭代次數，則終止迭代；否則令構造新的等式約束優化問題，轉3）。

優化迭代后，求解結果即為術前規劃的最優解，橈骨遠端畸形骨段在最優截骨位置繞復位旋轉軸進行旋轉復位，旋轉角度為優化計算的最優復位角。

2 實驗與結果分析

2.1 實驗數據集及環境

實驗數據為6 例一側畸形、對側健康的患者橈骨CT，來源于武漢大學中南醫院。圖像層內間距0.78 mm，層間間距1.25 mm，圖像強度單位為HU，對骨骼CT 進行重建處理重建后的三維模型格式為STL。

基于Windows 10 系統，采用Visual Studio 2019 開發平臺，配置VTK8.2.0，PCL1.11.0 及QT5.14.0，開發出三維截骨自動規劃原型系統，具有模型配準、切割及變換、骨畸形區域自動計算、橈骨成角畸形自動截骨矯形等功能。實驗計算機硬件配置為Intel Core i7-9700F CPU @3.70 GHz，16 GB內存。

2.2 復位精度評估

選取6 例畸形橈骨模型，如圖7 所示。算法驗證組設置內點算法的初始懲罰因子大小r0=1，懲罰因子的縮小系數c=0.5，精度ε=10-4，迭代次數為200，通過系統自動計算最優的截骨位置及最優的復位角。

圖7 畸形骨與對側鏡像健康骨原始模型Fig.7 Original model of deformed bone and contralateral mirror healthy bone

使用Miyake 等［6］提出的上肢畸形愈合骨折的三維截骨方法作為手動規劃實驗對照組。該方法對橈骨遠端和遠端關節整體配準對齊，計算畸形軸，醫師根據經驗在適當部位實施截骨術?；屋S沿凹側時，手動選擇兩個截骨平面切除楔形物，使遠端骨段繞畸形軸旋轉完成閉合楔形截骨術?；屋S沿凸側時，手動選擇截骨平面，切開楔形骨完成開放楔形截骨術。

使用文獻［11］中提出的楔形三維截骨方法作為自動規劃實驗對照組。該方法對前臂畸形骨與遠、近端復位模板進行配準，解耦配準后的變換矩陣，獲取畸形軸位置及矯形角度，以畸形軸位置作為截骨平面位置，將切割后的遠端骨段自動繞畸形軸旋轉所需角度，以實現矯正。

傳統上，畸形根據比較量化，可與健康的對側進行比較，因此，計算算法驗證組和實驗對照組中每一例畸形矯正后的遠端關節解剖區域配準精度誤差，如表2、3 所示。

表2 開放楔形關節解剖區域配準RMSETab.2 RMSE of anatomical area registration of open wedge joint

表3 閉合楔形關節解剖區域配準RMSETab.3 RMSE of anatomical area registration of closed wedge joint

圖8 是6 例畸形橈骨實例矯正效果對比，Case 1～Case 4為開放楔形截骨，Case 5～Case 6 為閉合楔形截骨。其中：A表示手工規劃組復位效果，B 為自動規劃對照組復位效果，C為本文算法的復位效果，D 為手工規劃組的截骨位置及矯正角度大小效果，E 為自動規劃組的截骨位置及矯正角度大小效果，F 為本文算法的截骨位置及矯正角度大小效果。

圖8 三維規劃實例圖Fig.8 Examples of 3D planning

由表2、3 可知，與手工規劃相比，每一例畸形骨的自動規劃誤差均小于手工規劃誤差，并且自動規劃的關節解剖區域配準RMSE 比手工規劃降低0.093 0～0.421 9 mm，平均降低0.211 4 mm，結果表明自動規劃的復位精度矯正效果均優于手動規劃。與自動規劃組相比，Case 1 中本文算法復位精度優于自動規劃組，而Case 2～Case 6 自動規劃組復位精度優于本文算法。

由圖8 分析可知，Case 2～Case 6 中，自動規劃組計算的旋轉軸遠離皮質骨邊緣，矯正后的骨頭無法明確定義楔形類型，部分重疊部分張開，如圖8 中E2～E6，這種情況可以在術前規劃時進行模擬，但在臨床上不可行，因為旋轉軸離皮質骨邊緣過遠，骨頭并未完全切斷，因此截骨后的遠端骨段無法繞剩余骨橋進行旋轉。本文算法計算的旋轉軸位置均處于畸形骨皮質骨表面輪廓上，在截骨復位的過程中避免了截骨處的嵌合，因此能在滿足臨床要求的前提下，明確定義開放楔形和閉合楔形截骨方式，并精準復位。如圖8 中F1～F4可直觀看出為開式楔形截骨，F5～F6 為閉式楔形截骨。Case 1 中，自動規劃組計算的復位旋轉軸位置離皮質骨邊緣較近，此時矯正的部分出現重疊范圍較小，如E1，此時可以明確定義為開放楔形并進行矯正復位，但是它的復位精度低于本文算法。

圖9 為復位旋轉軸位置示意圖。

圖9 復位旋轉軸位置示意圖Fig.9 Position diagram of reset rotation axis

綜上，根據實驗結果可知本文算法復位精度及可行性更高。

3 結語

針對橈骨成角畸形，本文提出了一種三維楔形截骨矯形自動規劃方法。選擇6 例成角畸形矯形案例，與醫生手動截骨矯形進行了對比分析，并通過橈骨遠端關節解剖區域配準RMSE 大小驗證。根據實驗數據結果可知：本文算法相較于手動規劃，復位精度更高；與自動規劃方法相比，本文算法可明確楔形類型，且臨床可行性更高。本文算法的局限在于，依賴于獲取患者整個前臂雙側的CT 掃描，可能導致額外的CT 輻射暴露；依賴于骨骼解剖，忽略了術前規劃中軟組織缺失的影響。對于未來的研究，首先考慮健康的對側無法參與規劃的情況，預測畸形骨未受影響部分的正常形狀；其次考慮術前過程中軟組織的影響從而實現更復雜的模擬規劃。