三維橈骨成角楔形截骨術(shù)前自動(dòng)規(guī)劃算法

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

（1.武漢理工大學(xué) 機(jī)電工程學(xué)院，武漢 430070；2.武漢大學(xué)中南醫(yī)院 骨科，武漢 430071）

0 引言

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

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

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

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

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

基于上述問題，針對(duì)橈骨成角畸形，本文提出一種橈骨楔形截骨矯形術(shù)前自動(dòng)規(guī)劃方法。首先，利用對(duì)側(cè)鏡像骨模型補(bǔ)償差異后的模型作為矯形復(fù)位模板，通過迭代最近點(diǎn)（Iterative Closest Point，ICP）算法對(duì)橈骨近端和遠(yuǎn)端關(guān)節(jié)進(jìn)行配準(zhǔn)，自動(dòng)計(jì)算畸形區(qū)域；其次，通過對(duì)橈骨遠(yuǎn)端關(guān)節(jié)的解剖區(qū)域進(jìn)行加權(quán)配準(zhǔn)，確定三維畸形軸的方向；隨后，調(diào)整畸形骨使它與世界坐標(biāo)系Y軸平行，并將畸形區(qū)域映射至XOZ平面，其中的畸形輪廓曲線采用三次樣條插值法確定，從而確定復(fù)位旋轉(zhuǎn)軸方位；最后，將關(guān)節(jié)解剖區(qū)域?qū)R作為截骨矯形的優(yōu)化目標(biāo)，以畸形區(qū)域?yàn)榻毓瞧矫婕s束、復(fù)位角的旋轉(zhuǎn)范圍為復(fù)位角約束，通過單目標(biāo)優(yōu)化算法進(jìn)行迭代求解，自動(dòng)生成最優(yōu)的橈骨楔形截骨矯形術(shù)前規(guī)劃方案。

1 楔形截骨矯形術(shù)前自動(dòng)規(guī)劃算法

1.1 算法流程

三維橈骨成角楔形截骨術(shù)前自動(dòng)規(guī)劃算法的總流程如圖1 所示，具體步驟如下：

圖1 算法流程Fig.1 Algorithm flow

1）將CT 數(shù)據(jù)進(jìn)行三維重建，生成畸形和對(duì)側(cè)健康橈、尺骨模型；

2）矢狀面作為對(duì)稱面，將對(duì)側(cè)健康的骨模型進(jìn)行鏡像，隨后將鏡像得到的健康模型用作畸形骨矯形復(fù)位的模板；

3）采用迭代最近點(diǎn)配準(zhǔn)算法對(duì)齊兩側(cè)橈骨和尺骨，計(jì)算兩側(cè)差異量并補(bǔ)償，重建新的矯形復(fù)位模板作為參考骨；

4）計(jì)算畸形橈骨的畸形區(qū)域范圍；

5）采用基于關(guān)節(jié)解剖區(qū)域權(quán)重的配準(zhǔn)方法，計(jì)算橈骨三維畸形軸方向；

6）畸形骨變換至與世界坐標(biāo)系Y軸平行，根據(jù)畸形區(qū)域范圍計(jì)算畸形輪廓在XOZ平面的曲線及復(fù)位旋轉(zhuǎn)軸方位；

7）畸形輪廓曲線上任取一點(diǎn)構(gòu)建復(fù)位旋轉(zhuǎn)軸，將畸形區(qū)域范圍和復(fù)位旋轉(zhuǎn)角度范圍作為控制變量約束，以橈骨遠(yuǎn)端解剖區(qū)域配準(zhǔn)均方根誤差（Root Mean Square Error，RMSE）最小為目標(biāo)設(shè)計(jì)目標(biāo)函數(shù)，構(gòu)建單目標(biāo)優(yōu)化模型；

8）優(yōu)化迭代計(jì)算出最優(yōu)的截骨平面和復(fù)位旋轉(zhuǎn)角度，自動(dòng)生成橈骨楔形截骨矯形術(shù)前規(guī)劃方案。

1.2 補(bǔ)償雙側(cè)差異

1）通過CT 數(shù)據(jù)創(chuàng)建畸形和健康側(cè)橈骨和尺骨的三維三角形表面模型。

2）將三維世界坐標(biāo)系中的XOZ平面作為對(duì)稱平面，獲取健康模型中每一個(gè)坐標(biāo)點(diǎn)關(guān)于XOZ平面的對(duì)稱點(diǎn)，根據(jù)點(diǎn)云數(shù)據(jù)生成對(duì)側(cè)健康前臂的鏡像模型。

3）采用雙側(cè)補(bǔ)償?shù)姆椒ǎ?8］糾正左右臂橈尺骨之間的長度差，提高術(shù)前規(guī)劃的復(fù)位精度。具體步驟如下：

a）計(jì)算橈骨和尺骨模型的質(zhì)心，并確定三條相互垂直的主軸線通過這兩個(gè)質(zhì)心。其中具有最小轉(zhuǎn)動(dòng)慣量的骨軸即為重力軸，對(duì)應(yīng)橈骨的主旋轉(zhuǎn)軸即橈骨長軸LR，對(duì)應(yīng)尺骨的主旋轉(zhuǎn)軸即尺骨長軸LU，如圖2 所示。

圖2 補(bǔ)償雙側(cè)差異Fig.2 Compensation for bilateral differences

b）橈骨近端和遠(yuǎn)端配準(zhǔn)對(duì)齊得到矩陣Mpr和Mdr，由這兩個(gè)矩陣可獲得修正矩陣使橈骨遠(yuǎn)端相對(duì)于近端重新定位。同時(shí)，尺骨近端和遠(yuǎn)端配準(zhǔn)對(duì)齊得到矩陣Mpu和Mdu，由這兩個(gè)矩陣可獲得補(bǔ)償矩陣

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

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

1.3 畸形區(qū)域

在優(yōu)化過程中計(jì)算截骨位置時(shí)，采用一種畸形區(qū)域的自動(dòng)計(jì)算方法［13］能夠明確畸形骨偏離參考骨的具體位置，縮小搜索空間，同時(shí)減少計(jì)算量和時(shí)間成本。

1）利用ICP 算法對(duì)畸形及參考骨進(jìn)行近、遠(yuǎn)端配準(zhǔn)，計(jì)算兩端的配準(zhǔn)點(diǎn)集；

2）沿畸形骨模型的骨長軸將配準(zhǔn)后的畸形骨模型和參考骨模型的點(diǎn)集以合適步長進(jìn)行空間離散處理；

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

其中：PSD表示盒內(nèi)畸形骨模型點(diǎn)集；PSN表示盒內(nèi)參考骨模型點(diǎn)集；|PSD|表示集合PSD的基數(shù)。自動(dòng)計(jì)算出的畸形函數(shù)圖像如圖3 所示。

圖3 骨畸形函數(shù)Fig.3 Bone deformity function

4）根據(jù)外科醫(yī)生建議，以偏離RMSE 最小值15%為畸形區(qū)域，計(jì)算畸形骨的畸形區(qū)域如圖4 中所示。

圖4 骨畸形區(qū)域Fig.4 Area of bone deformity

1.4 三維畸形軸

1.4.1 基于權(quán)重的關(guān)節(jié)局部區(qū)域配準(zhǔn)

在臨床上，橈骨遠(yuǎn)端由于與腕關(guān)節(jié)及尺骨遠(yuǎn)端關(guān)節(jié)相連，它的關(guān)節(jié)面的對(duì)齊比長骨區(qū)域的對(duì)齊更重要。因此，采用基于權(quán)重的關(guān)節(jié)解剖區(qū)域配準(zhǔn)方法［13］對(duì)畸形骨與參考骨遠(yuǎn)端的標(biāo)記區(qū)域進(jìn)行精準(zhǔn)復(fù)位對(duì)齊，以減小粗略對(duì)齊產(chǎn)生的誤差。

1）在遠(yuǎn)端關(guān)節(jié)解剖區(qū)域選取7 個(gè)解剖點(diǎn)，并對(duì)7 個(gè)解剖點(diǎn)確定的解剖區(qū)域設(shè)置不同權(quán)重，權(quán)重取值范圍為［0，1］，具體分配情況見表1。

表1 解剖點(diǎn)序號(hào)對(duì)應(yīng)解剖位置名稱及權(quán)重值Tab.1 Anatomic point serial number and corresponding anatomic position name and weight value

2）采用KDtree 算法［22］在7 個(gè)解剖點(diǎn)鄰域內(nèi)搜索50 個(gè)點(diǎn)。ICP 配準(zhǔn)的評(píng)價(jià)指標(biāo)［23］：待配準(zhǔn)點(diǎn)集與參考點(diǎn)集之間經(jīng)過矩陣變換后，點(diǎn)集數(shù)據(jù)之間滿足某種度量準(zhǔn)則下的最優(yōu)匹配，即平均距離差值最小。因此，提出一種基于區(qū)域權(quán)重的點(diǎn)云配準(zhǔn)誤差評(píng)價(jià)指標(biāo)如下：其中：K表示橈骨遠(yuǎn)端解剖區(qū)域的總數(shù)；N=50 為最鄰近點(diǎn)對(duì)個(gè)數(shù)；wi表示各個(gè)解剖區(qū)域的權(quán)重值；R表示從畸形骨解剖點(diǎn)集變換至參考骨解剖點(diǎn)集的旋轉(zhuǎn)矩陣；pij表示畸形骨解剖區(qū)域i中的第j個(gè)最近點(diǎn)，qij表示參考骨鏡像模型解剖區(qū)域i中的第j個(gè)最近點(diǎn)。

1.4.2 三維畸形軸方向

1）畸形骨和參考骨遠(yuǎn)端通過關(guān)節(jié)整體區(qū)域配準(zhǔn)，獲得變換矩陣Md；畸形骨和參考骨近端通過基于權(quán)重的關(guān)節(jié)局部區(qū)域配準(zhǔn)，獲得變換矩陣Mp，由兩個(gè)變換矩陣產(chǎn)生校正矩陣

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

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

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

1.5 畸形輪廓曲線和復(fù)位旋轉(zhuǎn)軸方位

1.5.1 旋轉(zhuǎn)軸方向向量

確定畸形軸的方向，構(gòu)建同向的單位向量并進(jìn)行變換，具體步驟如下：

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

2）計(jì)算方向向量V變換到與世界坐標(biāo)系Y軸平行的變換矩陣MY；

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

1.5.2 畸形輪廓曲線

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

在每一個(gè)小區(qū)間(zi-1，zi)(i=1，2，…，n)內(nèi)，S(z)分別為三次多項(xiàng)式函數(shù)：

將數(shù)據(jù)節(jié)點(diǎn)和指定的首尾端點(diǎn)條件代入矩陣方程求解可得樣條曲線系數(shù)ai、bi、ci、di，即可計(jì)算出曲線方程S1(z)，曲線方程S2(z)。畸形輪廓曲線圖如圖6 所示。

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

1.5.3 復(fù)位旋轉(zhuǎn)軸方位

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

1.6 優(yōu)化方法

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

1.6.1 變量約束

變量約束包括截骨位置約束和復(fù)位旋轉(zhuǎn)軸旋轉(zhuǎn)的角度約束。

截骨位置約束是指在截骨過程中截骨平面可變化的范圍。截骨平面決定截骨的大小和位置以及截骨的方向，不同的截骨位置會(huì)導(dǎo)致復(fù)位誤差不同。優(yōu)化開始前，通過骨畸形自動(dòng)診斷確定截骨平面大致范圍，即為截骨平面位置約束。

復(fù)位旋轉(zhuǎn)軸旋轉(zhuǎn)的角度約束是遠(yuǎn)端骨段繞復(fù)位旋轉(zhuǎn)軸旋轉(zhuǎn)至與參考骨模型遠(yuǎn)端的解剖區(qū)域?qū)R，即為復(fù)位旋轉(zhuǎn)軸旋轉(zhuǎn)的角度約束。旋轉(zhuǎn)的角度θ變化范圍為(-45°，45°)。

1.6.2 目標(biāo)函數(shù)

橈骨成角畸形截骨術(shù)的術(shù)前優(yōu)化目標(biāo)為橈骨遠(yuǎn)端的關(guān)節(jié)局部區(qū)域配準(zhǔn)RMSE，用于測量配準(zhǔn)精度，精細(xì)控制復(fù)位對(duì)齊。具體描述如下：

橈骨遠(yuǎn)端關(guān)節(jié)解剖區(qū)域的配準(zhǔn)精度：

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

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

1.6.3 內(nèi)點(diǎn)算法

將實(shí)際問題轉(zhuǎn)化為算法求解的形式，具體步驟如下：

1）建立數(shù)學(xué)模型。式（13）為目標(biāo)函數(shù)，式（14）為不等式約束集合：

2）引入松弛因子，構(gòu)造對(duì)數(shù)障礙函數(shù)作為懲罰函數(shù)，生成等價(jià)于原模型的優(yōu)化問題，并將不等式約束問題轉(zhuǎn)化為等式約束問題。

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

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

4）構(gòu)造拉格朗日函數(shù)，處理等式約束。

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

5）利用迭代法求解庫恩塔克條件（式（17）），獲得最優(yōu)點(diǎn)。

求解過程中，當(dāng)拉格朗日函數(shù)的Hessian 矩陣正定時(shí)，目標(biāo)函數(shù)在迭代點(diǎn)附近鄰域內(nèi)是嚴(yán)格凸函數(shù)，算法采用牛頓法。當(dāng)Hessian 矩陣半正定時(shí)，目標(biāo)函數(shù)在迭代點(diǎn)附近鄰域內(nèi)是凹函數(shù)，牛頓法失效，算法轉(zhuǎn)為基于信賴域的共軛梯度法求解。逼近問題（15）在迭代附近不是局部凸時(shí)，直接采用基于信賴域的共軛梯度法。

6）若‖x?(r(k))-x?(r(k-1)) ‖≤ε或達(dá)到迭代次數(shù)，則終止迭代；否則令構(gòu)造新的等式約束優(yōu)化問題，轉(zhuǎn)3）。

優(yōu)化迭代后，求解結(jié)果即為術(shù)前規(guī)劃的最優(yōu)解，橈骨遠(yuǎn)端畸形骨段在最優(yōu)截骨位置繞復(fù)位旋轉(zhuǎn)軸進(jìn)行旋轉(zhuǎn)復(fù)位，旋轉(zhuǎn)角度為優(yōu)化計(jì)算的最優(yōu)復(fù)位角。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 實(shí)驗(yàn)數(shù)據(jù)集及環(huán)境

實(shí)驗(yàn)數(shù)據(jù)為6 例一側(cè)畸形、對(duì)側(cè)健康的患者橈骨CT，來源于武漢大學(xué)中南醫(yī)院。圖像層內(nèi)間距0.78 mm，層間間距1.25 mm，圖像強(qiáng)度單位為HU，對(duì)骨骼CT 進(jìn)行重建處理重建后的三維模型格式為STL。

基于Windows 10 系統(tǒng)，采用Visual Studio 2019 開發(fā)平臺(tái)，配置VTK8.2.0，PCL1.11.0 及QT5.14.0，開發(fā)出三維截骨自動(dòng)規(guī)劃原型系統(tǒng)，具有模型配準(zhǔn)、切割及變換、骨畸形區(qū)域自動(dòng)計(jì)算、橈骨成角畸形自動(dòng)截骨矯形等功能。實(shí)驗(yàn)計(jì)算機(jī)硬件配置為Intel Core i7-9700F CPU @3.70 GHz，16 GB內(nèi)存。

2.2 復(fù)位精度評(píng)估

選取6 例畸形橈骨模型，如圖7 所示。算法驗(yàn)證組設(shè)置內(nèi)點(diǎn)算法的初始懲罰因子大小r0=1，懲罰因子的縮小系數(shù)c=0.5，精度ε=10-4，迭代次數(shù)為200，通過系統(tǒng)自動(dòng)計(jì)算最優(yōu)的截骨位置及最優(yōu)的復(fù)位角。

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

使用Miyake 等［6］提出的上肢畸形愈合骨折的三維截骨方法作為手動(dòng)規(guī)劃實(shí)驗(yàn)對(duì)照組。該方法對(duì)橈骨遠(yuǎn)端和遠(yuǎn)端關(guān)節(jié)整體配準(zhǔn)對(duì)齊，計(jì)算畸形軸，醫(yī)師根據(jù)經(jīng)驗(yàn)在適當(dāng)部位實(shí)施截骨術(shù)?；屋S沿凹側(cè)時(shí)，手動(dòng)選擇兩個(gè)截骨平面切除楔形物，使遠(yuǎn)端骨段繞畸形軸旋轉(zhuǎn)完成閉合楔形截骨術(shù)?；屋S沿凸側(cè)時(shí)，手動(dòng)選擇截骨平面，切開楔形骨完成開放楔形截骨術(shù)。

使用文獻(xiàn)［11］中提出的楔形三維截骨方法作為自動(dòng)規(guī)劃實(shí)驗(yàn)對(duì)照組。該方法對(duì)前臂畸形骨與遠(yuǎn)、近端復(fù)位模板進(jìn)行配準(zhǔn)，解耦配準(zhǔn)后的變換矩陣，獲取畸形軸位置及矯形角度，以畸形軸位置作為截骨平面位置，將切割后的遠(yuǎn)端骨段自動(dòng)繞畸形軸旋轉(zhuǎn)所需角度，以實(shí)現(xiàn)矯正。

傳統(tǒng)上，畸形根據(jù)比較量化，可與健康的對(duì)側(cè)進(jìn)行比較，因此，計(jì)算算法驗(yàn)證組和實(shí)驗(yàn)對(duì)照組中每一例畸形矯正后的遠(yuǎn)端關(guān)節(jié)解剖區(qū)域配準(zhǔn)精度誤差，如表2、3 所示。

表2 開放楔形關(guān)節(jié)解剖區(qū)域配準(zhǔn)RMSETab.2 RMSE of anatomical area registration of open wedge joint

表3 閉合楔形關(guān)節(jié)解剖區(qū)域配準(zhǔn)RMSETab.3 RMSE of anatomical area registration of closed wedge joint

圖8 是6 例畸形橈骨實(shí)例矯正效果對(duì)比，Case 1～Case 4為開放楔形截骨，Case 5～Case 6 為閉合楔形截骨。其中：A表示手工規(guī)劃組復(fù)位效果，B 為自動(dòng)規(guī)劃對(duì)照組復(fù)位效果，C為本文算法的復(fù)位效果，D 為手工規(guī)劃組的截骨位置及矯正角度大小效果，E 為自動(dòng)規(guī)劃組的截骨位置及矯正角度大小效果，F(xiàn) 為本文算法的截骨位置及矯正角度大小效果。

圖8 三維規(guī)劃實(shí)例圖Fig.8 Examples of 3D planning

由表2、3 可知，與手工規(guī)劃相比，每一例畸形骨的自動(dòng)規(guī)劃誤差均小于手工規(guī)劃誤差，并且自動(dòng)規(guī)劃的關(guān)節(jié)解剖區(qū)域配準(zhǔn)RMSE 比手工規(guī)劃降低0.093 0～0.421 9 mm，平均降低0.211 4 mm，結(jié)果表明自動(dòng)規(guī)劃的復(fù)位精度矯正效果均優(yōu)于手動(dòng)規(guī)劃。與自動(dòng)規(guī)劃組相比，Case 1 中本文算法復(fù)位精度優(yōu)于自動(dòng)規(guī)劃組，而Case 2～Case 6 自動(dòng)規(guī)劃組復(fù)位精度優(yōu)于本文算法。

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

圖9 為復(fù)位旋轉(zhuǎn)軸位置示意圖。

圖9 復(fù)位旋轉(zhuǎn)軸位置示意圖Fig.9 Position diagram of reset rotation axis

綜上，根據(jù)實(shí)驗(yàn)結(jié)果可知本文算法復(fù)位精度及可行性更高。

3 結(jié)語

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