朱新成,何坤金,3*,倪 娜,郝 博
(1.河海大學物聯網工程學院,江蘇常州 213022;2.常州市圖形圖像與骨科植入物數字化技術重點實驗室,江蘇常州 213022;3.疏浚技術教育部工程研究中心,江蘇常州 213022)
隨著我國人口老齡化的加劇和生活中各類交通事故的頻發,骨折患者不斷增加。目前,外科手術中最常見的骨科植入物有接骨板、髓內釘、鋼針等,其中接骨板的設計是植入物設計的核心[1]。接骨板作為臨床上應用最為普遍的內植入物之一,其作用是保持骨折端的復位狀態,提供良好的穩定性[2],因此接骨板在手術中貼合位置的準確性對患者的術后恢復起著關鍵作用[3]。由于術前醫生只能大致判斷接骨板的貼合位置,并在術中需要對接骨板位置進行多次調整,這無疑增加了手術的難度,延長了手術的時間。因此,提前計算出接骨板的最佳貼合位置能夠指導醫生更好地進行術前規劃,這對于計算機輔助醫生進行高效診療有著重要的意義。
近年來,數字骨科[4]不斷發展,成為推動骨科臨床發展的重要動力。Schmutz 等[5]通過計算機在全虛擬環境下對脛骨遠端內側接骨板的貼合性進行評估,并提出了在不同部位接骨板到骨表面之間間隙的最低要求;吳旭華等[6]利用接骨板與骨骼表面之間間隙的體積和接骨板內曲面面積之間的關系來對不同接骨板的貼合性進行比較;Harith等[7]提出使用多個貼合指標以及接骨板與骨面的最大距離衡量臨床治療結果,并開發了一種用于自動鋼板定位和貼合性評估的軟件工具,以達到植入物設計驗證和優化的目的;Malekani 等[8]分析了316L和Ti-6Al-4V不銹鋼接骨板并進行壓力性能測試,為手術中貼合性不佳且需要手動彎曲調整的接骨板提供了一個指導性標準。貼合性良好的接骨板可以更好地幫助醫生借助接骨板對骨折進行復位,減少接骨板對軟組織的影響[9]。上述研究提出了接骨板貼合指標,衡量了部分廠家的接骨板貼合性,但并未提供有效的方式指導醫生確定接骨板貼合的最佳位置。針對接骨板的研究有了一定的成果[10-12],提出了個性化和系列化設計,為未來接骨板的發展趨勢提供了新思路。由于技術水平和成本投入的限制,目前國內接骨板生產廠家仍是批量生產,這就使得即使接骨板的型號分得足夠細,醫生根據經驗選取完合適的接骨板在手術中仍然需要多次比對找到與骨骼表面貼合最佳的位置,至今針對接骨板貼合性的計算方法研究較少。針對點云配準技術,Besl 等[13]提出了迭代最近點(Iterative Closest Point,ICP)算法,由于操作簡單以及良好的配準效果得到了廣泛的應用,但是ICP 算法對于點云的初始位置要求很高且配準效率低;Mitra 等[14]提出了一種基于二次距離框架的配準算法,同時利用點與點和點與面之間的距離來達到較好的配準效果,解決了兩組點云初始距離遠的問題,但是該算法需要計算每個點的曲率,導致計算量大,配準時間長;王欣等[15]提取點云模型邊界的特征點,在初始配準的基礎上,運用改進的ICP 算法進行精確配準,雖然算法效率提高,但是精確度不高;陳學偉等[16]基于采樣一致性初始配準(SAmple Consensus Initial Aligment,SAC-IA)算法進行初始配準,再利用改進的ICP 算法精確配準,達到了較好的配準效果,但是配準時間較長;Liu 等[17]將K 均值聚類算法與改進權值比的ICP算法相結合,較傳統的ICP算法雖然降低了配準時間,但配準效果不佳。。
由于脛骨周圍的肌肉組織較少,骨頭迎面區域在受到旋轉用力或直接暴力時更容易發生骨折,因此,本文針對上述問題,以脛骨為例進行研究,在對斷骨復位的研究成果[18]基礎上,提出了一種基于改進ICP 算法的接骨板貼合性快捷計算方法。首先,選取脛骨表面的貼合區域,提取接骨板內曲面點云;然后,對兩組點云平滑、采樣處理,利用點云之間的特征關系初始配準;最后,提取接骨板特征關鍵點,采用K-維樹(KDimensional Tree,KD-Tree)搜索鄰近點對接骨板點云模型和脛骨區域點云模型執行ICP精確配準。
為了指導醫生能快捷有效地在斷骨上確定接骨板最佳貼合位置,本文提出了一種基于改進ICP 算法的接骨板貼合性快捷計算方法,并以脛骨為例進行實驗,算法流程如圖1所示。

圖1 本文方法流程Fig.1 Flowchart of the proposed method
根據設計流程,基于改進ICP 算法的接骨板貼合性快捷計算方法的主要步驟如下:
1)由醫生指導通過鼠標框選出貼合區域,再計算接骨板點云的法向量,利用內、外曲面法向量之間夾角大于90°的特征提取內曲面點云。
2)對兩組點云模型平滑處理去除噪點,再對兩組點云模型格點采樣以簡化模型,利用點云之間的特征關系進行初始配準。
3)提取接骨板內曲面點云模型的邊界及內部的特征關鍵點,采用KD-Tree 搜索鄰近點對接骨板點云模型和脛骨區域點云模型執行ICP精確配準。
本文通過接骨板內曲面與脛骨表面距離來衡量接骨板貼合性,分別對復位后的脛骨和接骨板模型進行預處理,其中包括醫生指導在脛骨表面選取貼合區域和利用接骨板表面法向量之間的夾角對接骨板內曲面提取。
由于外科手術前醫生只能通過斷骨CT 圖像大致確定接骨板的貼合位置,根據文獻[18-19]已有的斷骨復位方法,對復位后的脛骨與接骨板進行匹配,醫生可以在術前規劃時以交互的方式在脛骨表面選取大致的貼合區域,選取的區域如圖2(a)所示。

圖2 脛骨和接骨板模型Fig.2 Model of tibia and orthopedic plate
接骨板模型由生產廠家掃描真實接骨板得到,由內曲面、外曲面和側面組成,為了提高接骨板與脛骨表面匹配度,接骨板在脛骨上的貼合性用接骨板內曲面與脛骨表面之間的距離衡量,因此,需對接骨板點云模型的內曲面進行提取。根據接骨板內、外曲面和側面中包含點的法向量差異較大的特點,可利用點的法向量之間的夾角來獲取接骨板內曲面,算法流程如下:
1)接骨板內曲面任意一點p0,計算出該點的法向量V0;
2)計算接骨板點云模型點(p1,p2,…,pn)的法向量(V1,V2,…,Vn)與點p0的法向量V0的夾角(a1,a2,…,an);
3)若點pi的法向量Vi與V0的夾角ai小于閾值T,則將點pi記為接骨板內曲面的點。
接骨板模型如圖2(b)所示,提取出的接骨板內曲面點云如圖2(c)所示。
在對接骨板點云模型進行掃描時,可能出現掃描儀光線對過于彎曲的表面無法完全覆蓋的情況,因此點云模型通常會出現點云不平滑甚至是空洞的狀況。為了能夠更加準確地表現點云特征,提高點云法線的準確性,采用移動最小二乘法(Moving Least Squares,MLS)[20]對周圍數據點進行高階多項式插值來平滑接骨板內曲面和脛骨選取區域點云模型。
在點云模型的一個局域子域U的擬合函數可以表示為:

式中:a(x)={a1(x),a2(x),…,an(x)}是坐標x的函數,為待求系數;p(x)={p1(x),p2(x),…,pn(x)}是一個階完備的多項式,通常表示為[1,u,v,u2,v2,uv]T。
加權離散L2范式[21]為:

式中;n表示影響區域節點的數目,f(x)表示擬合函數,yi是x=xi處的節點值,w(x-xi)是節點xi的權函數。為實現點云平滑處理,權函數通常選擇立方樣條權函數:

式中:x=;hi為第i個節點的權函數支持域的值;β為引入的影響系數范圍,通常取1.2~2.5。
為了降低計算量,節約計算時間,在保證接骨板貼合性不降低的前提下需要對兩組點云模型進行精簡。目前常見的點云精簡方法是對點云模型進行采樣,采樣方法有均勻采樣、幾何采樣、隨機采樣、格點采樣等。為了保證精簡后模型配準效果的有效性,采樣點需要滿足分布均勻且能較好地體現幾何特征的條件,因此對接骨板內曲面點云模型分別進行均勻采樣和格點采樣,采樣對比結果如圖3所示。

圖3 采樣結果對比Fig.3 Comparison of sampling results
接骨板內曲面點云個數為6 058,精簡后的個數為833,均勻采樣和格點采樣所耗費的時間分別為0.108 s 和0.049 s;從圖3 中可以看出,接骨板內曲面模型經過格點采樣后的點云分布較均勻采樣后的點云分布更為均勻且幾何特征表現得更加明顯。在本文實驗中,格點采樣的效率更高且采樣結果更加符合本文實驗的要求,所以本文采用格點采樣的方法來精簡點云模型。
計算點云模型各個頂點在x、y、z三個坐標軸上投影的最小值(lx,ly,lz)和最大值(mx,my,mz),過(lx,ly,lz)和(mx,my,mz)點,分別沿著x、y、z三個坐標軸方向構建一個長方體包圍盒,以一定的間距R分別沿著x、y、z三個方向對最小包圍盒進行等間距劃分,最小包圍盒內包含了若干個邊長為R的小方格,每個小方格內有若干個點,取距離小方格中心最近的點即為采樣點。由于脛骨點云模型相較于接骨板內曲面點云模型點云更加密集,幾何特征也更加豐富,為了在精簡模型的同時能夠盡可能多地保留幾何特征,對脛骨選取區域精簡的采樣半徑R1應小于對接骨板內曲面精簡的采樣半徑R2,脛骨選取區域和接骨板內曲面模型的精簡結果分別如圖4(a)、圖4(b)所示。

圖4 點云模型精簡結果Fig.4 Simplification results of point cloud model
接骨板最佳貼合位置的尋找主要通過接骨板內曲面點云模型和脛骨選取區域點云模型配準來實現,為了實現高精度、高效率的配準,本文將改進的配準算法分為初始配準和精確配準,根據兩組點云模型的特征采用SAC-IA 算法配準,能夠避免精確配準時陷入局部最優的情況;提取接骨板內曲面特征關鍵點,采用KD-tree 搜索鄰近點對兩組點云模型執行ICP配準。
本文初始配準利用了接骨板內曲面和脛骨選取區域點云模型之間點云的特征,首先分別計算兩組點云的特征直方圖(Fast Point Feature Histogram,FPFH);然后使用SAC-IA 算法[22]對接骨板內曲面和脛骨選取區域的點云模型進行初始配準,初始配準算法流程如下:
1)在接骨板內曲面點云模型中選取n個點作為采樣點,為了滿足每個采樣點FPFH 特征之間的差異性,采樣點之間的距離應滿足大于閾值dmin。
2)以點云模型中任意一點pi為圓心,計算半徑為r范圍內點云的特征直方圖,然后在脛骨選取區域點云模型中搜索與接骨板內曲面點云模型中采樣點FPFH 特征相似的點,將該點作為接骨板內曲面點云模型在脛骨選取區域內對應的點。
3)計算接骨板內曲面點云模型和其在脛骨選取區域對應的點云模型之間的剛體變換矩陣,懲罰函數表示為Es,使用Huber函數來衡量配準誤差,ei表示每對點之間的誤差,te為定義的誤差范圍。

重復上述三個步驟,對點云模型內所有的點進行計算,當誤差取最小時獲得接骨板內曲面點云模型和其在脛骨選取區域對應的點云模型之間的變換矩陣即為最優變換矩陣,脛骨選取區域點云模型與接骨板內曲面點云模型初始位置如圖5(a)所示,經過初始配準后的位置如圖5(b)所示。

圖5 初始配準結果Fig.5 Initial registration results
初始配準是將距離較遠且在不同坐標系下的兩組點云模型配對到同一坐標系下,并為精確配準提供較好的初始位置。精確配準采用改進的迭代最近點(ICP)算法,在空間點云搜索上利用KD-Tree的方式搜索,來進行計算,流程如圖6所示。

圖6 改進ICP算法流程Fig.6 Flowchart of improved ICP algorithm
改進的ICP算法步驟如下:
1)對接骨板內曲面點云模型采用改進的內部形態描述算子(Insic Shape Signatures,ISS)算法提取特征關鍵點。
2)采用KD-Tree 的方式在脛骨選取區域查找與接骨板內曲面點的最近點,得到兩組點云的初始對應關系。
3)根據步驟2)中的對應關系來計算兩組點云之間的變換矩陣,使得目標函數(式(9))的值最小,并得到接骨板內曲面點在脛骨選取區域對應的點。
4)根據收斂條件判斷兩組點云之間的距離d是否小于設定的閾值μ:若d大于μ,且未達到設置的迭代次數,返回步驟2)繼續迭代;否則迭代結束。
3.2.1 接骨板特征關鍵點提取
在進行精準匹配時,如果點云的規模龐大,則需要花費大量的時間進行計算,效率較低。為了在減少計算量的同時能保證計算的精確度,則需要對點云模型的特征關鍵點進行提取,本文針對接骨板內曲面點云模型進行特征關鍵點提取。
內部形態描述算子[23]能夠高度區分點云模型的特征向量,對三維模型的幾何特征進行編碼,完成高質量的點云配準。接骨板內曲面點云邊界作為曲面表達的重要幾何特征,ISS 算法提取特征關鍵點時并不能較好地保留點云邊界點。為此,在ISS 算法的基礎上進行改進,在保留點云邊界的基礎上,對接骨板內曲面點云模型進行特征點提取,接骨板內曲面特征關鍵點提取算法流程如下:
1)對接骨板內曲的點bi及其鄰域內的點(b1,b2,…,bn)做平面擬合,若向量bib1,bib2,…,bibn的方向偏向某一側且向量之間的夾角最大值θmax小于閾值θτ,則將該點記為邊界點。
2)計算接骨板內曲面點云模型中的每個點pi,以該點為中心、rd為半徑的范圍內點的權重為wi,wi與范圍內點的個數成反比,表示為:

3)計算rd范圍內所有的點pj(j=1,2,…,n)與點pi之間的協方差矩陣cov(pi):

4)計算每個點pi的協方差矩陣的特征值{λi,1,λi,2,λi,3}并按照降序排列;
5)設定閾值ε1和ε2,ε1和ε2的范圍在0~1,滿足λi,2/λi,1≤ε1和λi,3/λi,2≤ε2的點即為特征點。
6)重復上述步驟,直至完成對接骨板內曲面點云模型內所有點的計算,得到接骨板內曲面點云模型的特征關鍵點。
對精簡后的接骨板內曲面點云模型特征提取結果如圖7(b)所示,從圖中可以看出,提取的特征關鍵點能夠很好地描述接骨板內曲面的幾何特性,從而保證了接骨板內曲面與脛骨區域之間的高精度配準。

圖7 特征關鍵點提取Fig.7 Extraction of feature key points
3.2.2 ICP配準
迭代最近點算法的本質是基于最小二乘法的最優匹配算法[24],對于點云P中的n個點,使用一個3×3 旋轉矩陣Ro和一個三維向量t來描述剛體變換。對于點云P中的任意一點pi(xi,yi,zi)和經過變換后的點,兩個點滿足以下條件:

重復選擇對應的點qi來計算最優剛體變換矩陣T,直到滿足最佳匹配的收斂精度,即使目標函數Ep最小化。

在解決ICP 算法對初始位置要求高和點云過于龐大的問題后,ICP 算法雖然能得到較好的配準效果,但是該算法本身的計算效率不高。而KD-Tree[25]是一種分割K維數據空間的數據結構,可應用于多維空間關鍵數據的搜索,且搜索具有搜索速度快的優點,因此本文采用KD-Tree 的方式對目標點云的空間進行搜索。
對接骨板內曲面模型的特征關鍵點與脛骨選取區域進行精確配準,配準結果如圖8(b)所示,獲取最終的變換矩陣Rs和三維向量ts,對接骨板模型Plate作變剛體換得到變換后的接骨板模型Plate':

最終輸出結果如圖8(c)所示。

圖8 接骨板配準結果Fig.8 Registration results of orthopedic plate
對本文所提算法在Microsoft Visual Studio 2019 平臺進行實驗,系統運行所需軟硬件環境如下:操作系統Windows 7 及以上,內存2 GB或以上,CPU 2.1 GHz或以上。實驗數據來源于課題組對斷骨復位的研究和廠家掃描的接骨板模型數據,根據Schmutz等[5]提出的接骨板貼合性衡量標準,脛骨內側接骨板和脛骨之間的間隙至少滿足以下3 點要求:1)接骨板遠近端近關節處距離骨面不超過2 mm;2)接骨板中部距離骨面不超過6 mm;3)接骨板近端距離骨面不超過4 mm。
文獻[16-17]中都對傳統的ICP 算法進行了改進,為了驗證所提算法的實用性和有效性,將所提算法與ICP 算法、文獻[16]算法、文獻[17]算法進行比較,分別計算了接骨板遠近端近關節處最大距離S1、接骨板中部與骨面之間最大距離S2,以及接骨板近端與骨面最大距離S3。以此為基礎,在本文實驗中接骨板貼合性定義為:1)當接骨板平均間隙小于2 mm 時,對應的貼合指數大于0.5,接骨板的貼合性為優;2)當接骨板平均間隙大于3 mm時,接骨板貼合指數小于0.3,接骨板貼合性為差;3)貼合指數在0.4~0.5的接骨板貼合性為良,在0.3~0.4 的貼合性為中。對匹配時間、平均間隙、貼合指數和貼合性進行記錄,如表1所示。

表1 不同算法所用時間和貼合性對比Tab.1 Time cost and fit comparison of different algorithms
由表1 中可以看出,本文算法與所對比的三種算法所得結果都滿足Schmutz 等[5]提出的脛骨內側接骨板和脛骨之間的間隙三點要求,但傳統ICP 算法雖然能夠在較短的時間內得到結果,但是貼合性較低;文獻[16]算法雖然能得到較好的貼合性,但是計算時間過長;文獻[17]算法雖然能夠在短的時間內得到較好的結果,但與本文算法相比,本文算法所用時間更短且貼合性更佳。不同算法得到的最終輸出結果如圖9所示。

圖9 不同算法配準結果對比Fig.9 Registration results comparison of different algorithms
為了驗證所提算法的通用性,除了對脛骨近端受損和T形接骨板進行了實驗,還對脛骨中端和遠端受損分別選取兩種不同類型的接骨板進行驗證,實驗所需配準時間、平均間隙和接骨板貼合性如表2所示。

表2 不同受損類型和接骨板類型所用時間和貼合性對比Tab.2 Time cost and fit comparison of different types of damage and orthopedic plate
從表2 中可以看出,接骨板配準時間和貼合性主要與接骨板內曲面結構和接骨板選擇的型號相關,對于內曲面結構簡單的接骨板,如T 形、弧形等接骨板而言,所需的配準時間較短,接骨板的貼合性也較好;對于內曲面結構復雜的接骨板,如異形接骨板,所需的配準時間較長,平均間隙相對較大,但是仍在允許的最優間隙范圍內。不同類型的配準效果對比如圖10 所示,受損位置由矩形框標出,配準結果表明所選接骨板能夠在選區的區域準確地找到貼合性最佳的位置。

圖10 不同受損類型配準結果Fig.10 Registration results of different damage types
本文提出了一種基于改進ICP 算法的接骨板貼合性快捷計算方法,主要包括:選取貼合區域,提取接骨板內曲面點云;對兩組點云平滑、采樣,利用點云特征初始配準;提取接骨板特征關鍵點進行精確配準。本文方法的主要工作有:1)將配準過程分為初始配準和精確配準兩步,配準效果更佳;2)對ISS 算法改進,使其能夠同時獲取點云邊界和內部特征點;3)采用KD-tree 的方式搜索鄰近點以提高ICP 算法的計算速度。通過對三種不同受損類型的脛骨和接骨板進行配準,實驗結果表明本文所提算法能夠高效、準確地計算出接骨板在脛骨上的最佳貼合位置。通過與其他配準算法作對比可以得出,所提算法在保證匹配精確度的同時能夠有效縮短配準時間。本文方法為有效判斷接骨板貼合位置提供了新思路,能夠輔助醫生更直觀、高效地進行術前規劃。未來將建立接骨板樣本庫,根據醫生所選區域實現對接骨板型號自動選擇,以更加高效、準確地計算接骨板貼合位置。