蔣俊鋒,孫曉莉,鄧子越,黃 瑞,何坤金,陳正鳴,劉宏偉,張文璽

(1. 河海大學(xué) 物聯(lián)網(wǎng)工程學(xué)院,江蘇 常州 213022; 2. 河海大學(xué) 計(jì)算機(jī)與信息學(xué)院,南京 210024; 3. 南京醫(yī)科大學(xué)附屬常州第二人民醫(yī)院,江蘇 常州 213003; 4. 溧陽市人民醫(yī)院,江蘇 溧陽 213300)

0 引 言

面向骨折修復(fù)的計(jì)算機(jī)輔助術(shù)前規(guī)劃中,需要先分割三維碎骨再修復(fù)碎骨[1-2],如圖1所示. 這種處理方法通常費(fèi)時(shí)費(fèi)力,原因如下: 1) 骨折形態(tài)隨機(jī),碎骨會發(fā)生相互接觸；2) 斷裂邊界的曲率有時(shí)不顯著；3) 碎骨拼接時(shí),對于細(xì)小碎骨通?？梢詠G棄；4) 骨骼不僅包含相對光滑的外表面,還包含形態(tài)復(fù)雜的內(nèi)部松質(zhì)骨. 目前已有的軟件自動化程度較低,相關(guān)方法不能較好地支持這種隨機(jī)斷裂的碎骨分割和拼接.

本文考慮對輸入的碎骨三角網(wǎng)格模型,提取碎骨外表面,并在骨骼斷裂和碎骨塊接觸邊界對碎骨進(jìn)行分割,然后以每個(gè)碎骨外表面為單元進(jìn)行碎骨拼接,完成骨折高效修復(fù),從而有效提高骨折修復(fù)術(shù)前規(guī)劃的自動化程度和效率. 研究表明,碎骨分割與拼接兩者相輔相成,將對側(cè)健康骨模型作為先驗(yàn)?zāi)０?能有效提高分割和拼接的正確率. 基于此,本文提出一種基于模板與數(shù)據(jù)驅(qū)動的三維碎骨分割與拼接方法,通過模板匹配、特征深度學(xué)習(xí)等方法,實(shí)現(xiàn)分割-拼接迭代求精.

1 相關(guān)工作

碎骨修復(fù)主要涉及碎骨形狀分割和拼接兩方面內(nèi)容. 碎骨分割包含提取骨骼外表面及將相互接觸的碎骨分隔開; 碎骨分割后,需要通過剛體變換對碎骨調(diào)整位置進(jìn)行拼接.

1.1 形狀分割

三維模型分割的目標(biāo)是將模型分割成互不相交的部分. 若需對各分割部分進(jìn)行語義標(biāo)記,則這種分割稱為語義分割. 經(jīng)典三維分割方法來源于二維圖像分割,主要包括區(qū)域增長、分水嶺、迭代聚類、層次聚類和邊緣分割等方法[3]. 區(qū)域增長法是一種交互式算法,效率較高,但存在種子點(diǎn)和增長停止條件不易選擇的問題; 分水嶺算法存在分割邊界鋸齒形狀、過分割等問題,通常需要邊界光滑和分割區(qū)域合并等后處理操作; 迭代聚類需預(yù)先指定分割數(shù)量,對于三維模型的初始分割較敏感; 層次聚類是一種多分辨率的聚類分割方法; 基于邊緣分割方法試圖找到分割邊緣,但邊緣常不易封閉.

近年來,基于數(shù)據(jù)驅(qū)動的三維形狀分割研究備受關(guān)注. 根據(jù)有無訓(xùn)練數(shù)據(jù),可將分割方法分為有監(jiān)督學(xué)習(xí)[4-9]和無監(jiān)督學(xué)習(xí)[10-11]. 根據(jù)三維形狀表示方式,三維形狀分割可分為基于多視圖[4]、點(diǎn)云[5-6]、圖[7]、網(wǎng)格[8]和體素[9]等方法. 其中,多視圖方法利用傳統(tǒng)的二維圖像卷積神經(jīng)網(wǎng)絡(luò),基于點(diǎn)云、圖、網(wǎng)格、體素的神經(jīng)網(wǎng)絡(luò)方法分別研究對應(yīng)表示方式的卷積、池化等操作; 此外,還有一種方式是先手工提取網(wǎng)格模型特征,然后以網(wǎng)格模型特征為神經(jīng)網(wǎng)絡(luò)的輸入,通過卷積、池化操作進(jìn)一步提煉分割特征[12-13]. 現(xiàn)有數(shù)據(jù)驅(qū)動方法通常先利用神經(jīng)網(wǎng)絡(luò)訓(xùn)練(有監(jiān)督)或聚類(無監(jiān)督)進(jìn)行初步分割,再利用圖割法等傳統(tǒng)方法修正邊界,其中,利用神經(jīng)網(wǎng)絡(luò)或聚類的標(biāo)簽概率作為圖割法目標(biāo)函數(shù)的約束項(xiàng). Zhou等[14]先訓(xùn)練兩個(gè)神經(jīng)網(wǎng)絡(luò),用于初步分割網(wǎng)格及其標(biāo)記邊界,然后利用圖割法結(jié)合邊界引導(dǎo)進(jìn)一步實(shí)現(xiàn)有效分割；Xu等[13]提出了一種三維牙齒模型分割方法,先用訓(xùn)練的兩個(gè)神經(jīng)網(wǎng)絡(luò)分別對牙齒和牙齦以及牙齒之間作語義分割,然后利用圖割、主成分分析等方法修正分割結(jié)果；Shu等[10]利用一種無監(jiān)督學(xué)習(xí)的三維形狀分割方法,先通過自編碼器提取特征,然后對特征用高斯混合模型進(jìn)行分割,最后用圖割法修正分割.

三維骨骼模型的分割方法,通常利用模板對骨骼進(jìn)行分割[15-16]. 本文提出一種基于模板分割健康骨骼的方法[17]. 該方法在提取少量特征點(diǎn)的基礎(chǔ)上,先對骨骼網(wǎng)格模型進(jìn)行語義分割,然后對分割得到的模板,基于Laplace變形進(jìn)行非剛體配準(zhǔn),從而指導(dǎo)同類骨骼模型實(shí)現(xiàn)快速兼容性分割.

對于碎骨分割,Shadid等[18]提出了一種基于隨機(jī)分水嶺變換的方法實(shí)現(xiàn)了踝部碎骨CT圖像分割. 目前,面向三維碎骨分割的研究尚未見文獻(xiàn)報(bào)道. 經(jīng)CT三維重建后,三維骨骼模型內(nèi)部包含形態(tài)復(fù)雜的內(nèi)部組織，如圖1所示,碎骨相互接觸的邊界形態(tài)隨機(jī)復(fù)雜,同時(shí),分割過程還需去除細(xì)小碎骨. 因此,利用目前已有方法提取外表面并同時(shí)進(jìn)行精準(zhǔn)碎骨分割十分困難.

1.2 碎片拼接

碎骨復(fù)位屬于碎片拼接范疇,主要用于文物復(fù)原[19]和碎骨復(fù)位[20-25]. 按照碎片數(shù)量,碎片拼接可分為一對一拼接和多對多拼接兩類. 多對多拼接研究[23,26-27]著重于兩個(gè)以上碎片拼接,碎片全局拼接是一個(gè)帶約束的非線性組合優(yōu)化問題. 當(dāng)碎骨數(shù)量較多時(shí),搜索空間較大,搜索效率較低. 目前的全局拼接搜索方法包含貪心法[19]、分支限界法[28]等. 其中,基于圖的貪心搜索方法效果較好,但有時(shí)很難找到全局最優(yōu)解,尤其在缺失碎片數(shù)量較多時(shí),貪心搜索方法易產(chǎn)生錯(cuò)誤拼接. 一對一拼接研究[28-29]常用方法包括四點(diǎn)配準(zhǔn)法[30]、相關(guān)法[23]、隨機(jī)一致性采樣[22]等,其中四點(diǎn)配準(zhǔn)法效率較高. 按照斷裂邊界不同可分為基于斷裂線的復(fù)位[21,27]和基于斷裂面的復(fù)位[19,23]. 目前的方法多數(shù)依據(jù)斷裂面曲率變化特性識別斷裂面,以斷裂面或斷裂面的子區(qū)域?yàn)槠唇訂卧?對拼接單元進(jìn)行剛體配準(zhǔn). 按照是否利用模板,碎片拼接可分為有模板拼接[21-24]和無模板拼接[25,31-33]兩種. 對位對線的粗隆間骨折碎骨復(fù)位方法[33]屬于無模板拼接,該方法較好解決了包含骨骼軸線成角約束的骨骼拼接問題,為臨床粗隆間骨折復(fù)位術(shù)前規(guī)劃提高了自動化程度. 有模板拼接通常利用對側(cè)健康骨骼[21-24]或形狀統(tǒng)計(jì)模型[34]作為模板,其關(guān)鍵是構(gòu)建碎骨外表面與模板的對應(yīng)關(guān)系,尤其當(dāng)碎骨面積較小或者碎骨外表面特征不顯著時(shí),對應(yīng)關(guān)系計(jì)算較難實(shí)現(xiàn).

2 方法概述

本文方法總體流程如圖2所示,輸入為包含內(nèi)部組織且未分割的三維碎骨模型,以及該碎骨模型的對側(cè)健康完整骨骼,輸出為分割且拼接好的碎骨外表面集合. 總體采用分割-拼接迭代求精的思想實(shí)現(xiàn)碎骨分割和拼接.

圖2 方法總體流程Fig.2 Overall process of method

1) 碎骨分割: 由于人體骨骼外表面較相似,具有統(tǒng)計(jì)學(xué)特征,骨骼斷裂面和骨骼松質(zhì)骨形態(tài)較復(fù)雜,因此,本文考慮通過訓(xùn)練神經(jīng)網(wǎng)絡(luò)初步提取碎骨外表面; 然后通過碎骨與模板匹配后,用模板與碎骨表面之間距離差異以及碎骨頂點(diǎn)的曲率變化信息作為判據(jù),進(jìn)一步確定外表面邊界.

2) 碎骨拼接: 先將碎骨剛體配準(zhǔn)到模板,模板匹配關(guān)鍵問題在于構(gòu)建碎骨到模板之間的密集對應(yīng)關(guān)系,本文通過訓(xùn)練神經(jīng)網(wǎng)絡(luò)回歸學(xué)習(xí)骨骼的表面特征; 然后將所有初步提取的碎骨外表面整體作為神經(jīng)網(wǎng)絡(luò)輸入,預(yù)測每個(gè)碎片的表面特征,同時(shí)預(yù)測模板的表面特征,通過特征匹配構(gòu)建碎片與模板的對應(yīng)關(guān)系,完成碎骨碎片與模板對齊; 最后用ICP(iterative closest point)[35]方法修正對齊結(jié)果; 在每次迭代過程中需重新計(jì)算碎骨碎片特征,可通過標(biāo)簽傳遞將對側(cè)骨模板特征映射到碎骨碎片.

3) 分割與拼接迭代求精: 碎骨分割、碎骨表面特征計(jì)算、碎骨拼接這三者相互依賴、相輔相成. 拼接結(jié)果越好,碎片整體越接近于一個(gè)健康骨骼的形態(tài),使碎骨表面的特征回歸精準(zhǔn)度更高,促進(jìn)模板匹配的精度,從而進(jìn)一步提高分割和拼接精度. 因此本文考慮三者不斷迭代,逐步求精.

本文算法步驟如下:

輸入: 輸入模型為碎骨模型M和健康對側(cè)骨模型N(M為碎骨模型點(diǎn)云矩陣,N為對側(cè)骨模型點(diǎn)云矩陣),最大迭代次數(shù)P,收斂閾值δ;

輸出: 修復(fù)的碎骨模型{RiMi+ti},Mi?M,Mi∩Mj=?,其中Mi為第i個(gè)碎骨點(diǎn)云矩陣,Ri和ti分別為第i個(gè)碎骨旋轉(zhuǎn)變換矩陣和平移變換矩陣；

3) While (迭代次數(shù)k

{

基于模板特征和碎骨特征將初步分割的碎骨配準(zhǔn)到模板,然后用ICP算法修正匹配;

break;

k=k+1;

}.

3 外表面提取與分割

圖3 碎骨外表面提取Fig.3 Extraction of external surface of broken bones

提取區(qū)分骨骼的外表面是一個(gè)二分類問題: 一類是碎骨外表面; 另一類是碎骨其他區(qū)域(包含骨骼內(nèi)部組織和斷裂區(qū)域). 本文選擇PointNet++[5]分割神經(jīng)網(wǎng)絡(luò),通過深度學(xué)習(xí)提取碎骨外表面. 訓(xùn)練數(shù)據(jù)為200個(gè)碎骨模型,人工標(biāo)注每個(gè)頂點(diǎn)上這兩類面的標(biāo)簽. 預(yù)測目標(biāo)碎骨模型的外表面時(shí),輸出每個(gè)頂點(diǎn)屬于內(nèi)外表面的概率,網(wǎng)絡(luò)預(yù)測正確率約為95%. 碎骨無斷裂的平滑區(qū)域?qū)儆谕獗砻娓怕瘦^大,而靠近邊界處的頂點(diǎn)以及碎骨內(nèi)部屬于外表面的概率較小. 碎骨外表面提取結(jié)果如圖3所示,其中圖3(B)為碎骨點(diǎn)云模型通過PointNet++神經(jīng)網(wǎng)絡(luò)預(yù)測后內(nèi)外表面分割結(jié)果,黃色標(biāo)記點(diǎn)為預(yù)測得到的碎骨外表面頂點(diǎn),本文方法選用網(wǎng)絡(luò)預(yù)測的中間結(jié)果,即內(nèi)外表面概率作為方法的輸入. 因此,為保證各碎骨的外表面不相鄰,本文取外表面概率范圍在[0.95,1]內(nèi)的點(diǎn),然后根據(jù)連通性,初步分割構(gòu)成各碎骨的外表面,如圖3(C)所示.

對碎骨內(nèi)外表面概率進(jìn)行預(yù)測的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示,其中n=3. 網(wǎng)絡(luò)結(jié)構(gòu)主要包含3個(gè)階段: 1) 多尺度分組抽象層,采用兩種尺度進(jìn)行分組,通過對點(diǎn)云模型下采樣提取局部特征; 2) 特征傳播層,對點(diǎn)云進(jìn)行上采樣,通過分層插值法傳播特征; 3) 多層感知機(jī),輸出最終描述符. 分割神經(jīng)網(wǎng)絡(luò)的損失函數(shù)采用交叉熵函數(shù),定義為

(1)

圖4 PointNet++預(yù)測模型Fig.4 PointNet++ prediction model

4 模板匹配

模板匹配的目標(biāo)是將每塊碎骨匹配到對側(cè)骨模板的對應(yīng)區(qū)域. 本文先利用特征匹配初步配準(zhǔn)碎骨和模板,然后用ICP算法進(jìn)一步修正配準(zhǔn). 其中關(guān)鍵問題是計(jì)算碎骨特征. 骨骼特征分為對側(cè)骨模板特征和初始碎骨特征. 本文采用深度神經(jīng)網(wǎng)絡(luò)PointNet++[5]回歸預(yù)測. 碎骨配準(zhǔn)到模板后,在每次迭代過程中都需計(jì)算碎骨特征,本文通過標(biāo)簽傳遞將對側(cè)骨模板特征映射到碎骨.

4.1 基于形狀統(tǒng)計(jì)模板的訓(xùn)練數(shù)據(jù)生成

為降低人工標(biāo)注成本,本文利用一種基于可變形模板生成訓(xùn)練數(shù)據(jù),該方法通過模板的參數(shù)化變形生成大量帶特征的訓(xùn)練數(shù)據(jù). 步驟如下: 先利用主成分分析法[34,36]構(gòu)建參數(shù)化的健康骨骼形狀統(tǒng)計(jì)模板; 然后取其平均模板(參數(shù)皆為0)做Laplace嵌入[37],即計(jì)算三維形狀Laplace矩陣的特征值和特征向量,按特征值由小到大排序,取前30個(gè)特征值對應(yīng)的特征向量,每個(gè)特征向量的維數(shù)等于模板骨骼的頂點(diǎn)數(shù),從而生成每個(gè)頂點(diǎn)上的30維向量特征; 最后使模板參數(shù)變形的模型拓?fù)湟恢?將平均骨骼頂點(diǎn)上的向量特征映射到變形骨骼的對應(yīng)頂點(diǎn),形成訓(xùn)練數(shù)據(jù).

4.2 基于上下文形狀的碎骨特征預(yù)測

利用完整骨骼進(jìn)行訓(xùn)練,則該網(wǎng)絡(luò)對于完整骨骼預(yù)測特征較好. 因此,為提高特征預(yù)測準(zhǔn)確度,本文將所有碎骨外表面作為一個(gè)整體預(yù)測特征,而不用單個(gè)碎骨作為網(wǎng)絡(luò)輸入.

對外表面碎骨特征進(jìn)行預(yù)測的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示,其中n=30,其損失函數(shù)定義為

(2)

5 計(jì)算碎片邊界

圖5 碎骨模型上的曲率場和距離場Fig.5 Curvature field and distance field on broken bone model

區(qū)域增長停止條件同時(shí)滿足碎片邊界局部曲率條件及碎片邊界與模板間的距離條件：

(3)

(4)

其中：p0為當(dāng)前增長點(diǎn),pi為p0測地距離為3 mm的鄰域范圍內(nèi)第i個(gè)點(diǎn),m為p0鄰域范圍內(nèi)的點(diǎn)數(shù)；Cp0為p0的平均曲率,Cpi為pi的平均曲率；δc為曲率閾值,δc=0.008時(shí)實(shí)驗(yàn)結(jié)果最優(yōu);qNi為點(diǎn)pi在模板N上的歐氏距離最近點(diǎn);δd為距離閾值,δd=0.01 mm時(shí)實(shí)驗(yàn)結(jié)果最優(yōu); 距離函數(shù)d定義為

(5)

6 實(shí) 驗(yàn)

6.1 骨骼數(shù)據(jù)

本文實(shí)驗(yàn)數(shù)據(jù)為8例脛骨骨折雙側(cè)掃描臨床病例,模板為碎骨對側(cè)骨鏡像后數(shù)據(jù),均來自于南京醫(yī)科大學(xué)附屬常州第二人民醫(yī)院. 其中4例為脛骨骨干骨折,4例為脛骨遠(yuǎn)端骨折,男女比例為5∶3.

三維碎骨與模板為CT影像雙側(cè)掃描數(shù)據(jù),采用Mimics與3-Matic[38]軟件進(jìn)行三維重建與預(yù)處理,處理過程如下: 利用Mimics軟件實(shí)現(xiàn)碎骨與對側(cè)骨的三維重建; 利用3-Matic軟件中網(wǎng)格優(yōu)化對碎骨和對側(cè)骨模型進(jìn)行重網(wǎng)格,在不改變模型外形結(jié)構(gòu)的情形下適當(dāng)簡化網(wǎng)格模型; 提取對側(cè)骨的表面模型并進(jìn)行鏡像操作,得到碎骨模板數(shù)據(jù).

本文實(shí)驗(yàn)中真實(shí)碎骨數(shù)據(jù)共208例,200例數(shù)據(jù)為PointNet++內(nèi)外表面預(yù)測神經(jīng)網(wǎng)絡(luò)的訓(xùn)練數(shù)據(jù),8例為實(shí)驗(yàn)測試數(shù)據(jù),測試數(shù)據(jù)不包含于訓(xùn)練數(shù)據(jù). 由于人體骨骼外表面形態(tài)相似,具有一定的統(tǒng)計(jì)特性,骨骼內(nèi)部松質(zhì)骨分布雜亂無章,因此200例訓(xùn)練數(shù)據(jù)能較好地區(qū)分碎骨內(nèi)外表面.

6.2 實(shí)驗(yàn)過程及結(jié)果

6.2.1 實(shí)驗(yàn)過程

本文用一例復(fù)雜碎骨實(shí)例驗(yàn)證本文碎骨分割與拼接方法的有效性. 圖6(A)為三維重建后的復(fù)雜碎骨,圖6(B)為碎骨模板,模板是碎骨對側(cè)骨鏡像后的外表面模型.

本文方法碎骨分割與拼接過程如下:

1) 將碎骨點(diǎn)云模型輸入到PointNet++分割神經(jīng)網(wǎng)絡(luò)中,預(yù)測得到碎骨點(diǎn)云模型內(nèi)外表面頂點(diǎn)概率,提取外表面概率大于0.95的頂點(diǎn),得到碎片互不相連的碎骨初步分割結(jié)果,如圖6(C)所示;

2) 將碎骨初步分割點(diǎn)云模型轉(zhuǎn)為網(wǎng)格模型,根據(jù)連通性將其自動分割為互不相連的碎骨碎片;

3) 將碎骨初步分割點(diǎn)云模型與模板點(diǎn)云模型輸入到PointNet++回歸預(yù)測網(wǎng)絡(luò)中,預(yù)測得到碎骨初始特征與模板特征,分別如圖6(D)和圖6(E)所示;

4) 根據(jù)特征將碎骨初匹配到模板上,再根據(jù)ICP算法進(jìn)一步修正匹配,如圖6(F)所示;

5) 根據(jù)碎骨與模板間的距離差異與碎骨表面曲率信息進(jìn)一步分割碎片邊界;

6) 通過標(biāo)簽傳遞將模板特征傳遞到碎片上,更新碎片特征,根據(jù)特征再次進(jìn)行匹配、分割,此過程迭代3次后收斂,碎片分割與匹配精度逐步提高.

標(biāo)記每個(gè)碎片最終分割結(jié)果在模板上的對應(yīng)點(diǎn),當(dāng)對下一個(gè)碎片進(jìn)行分割時(shí)將不再匹配模板上已標(biāo)記的區(qū)域. 對碎骨所有碎片依次執(zhí)行上述步驟,圖7為碎骨所有碎片3次迭代處理過程,碎片紅色標(biāo)記表示碎片分割邊界,隨著迭代次數(shù)的增加,碎片邊界逐步優(yōu)化,碎骨拼接結(jié)果接近健康骨骼狀態(tài).

圖6 碎骨分割過程Fig.6 Segmentation process of broken bones

圖7 碎骨3次迭代分割拼接過程Fig.7 Three iterations segmentation and splicing process of broken bones

本文方法分割與拼接迭代求精,隨著迭代次數(shù)的增加,分割錯(cuò)誤率與拼接錯(cuò)誤率逐步下降. 圖8(A)為該例碎骨隨迭代次數(shù)增加分割錯(cuò)誤率的變化情況,采用分類指標(biāo)中錯(cuò)誤率計(jì)算方法,定義為

error_rate=(FP+FN)/(P+N)，

(6)

其中FP表示算法分割錯(cuò)誤三角面片數(shù)量,FN表示算法欠分割三角面片數(shù)量,P+N表示該模型三角面片總數(shù)量. 圖8(B)為拼接錯(cuò)誤率的變化情況,拼接錯(cuò)誤率從平移相對誤差以及x,y,z軸旋轉(zhuǎn)相對誤差評估,定義為

(7)

其中n表示碎骨包含的碎片個(gè)數(shù),Δi表示第i個(gè)碎片的絕對誤差,即手工測量值與算法測量值的絕對差值,Li表示手工測量值. 分割錯(cuò)誤率與拼接錯(cuò)誤率在3次迭代后收斂,為保證算法效率,本文方法迭代次數(shù)最終選擇為3次.

圖8 碎骨分割(A)與碎骨拼接(B)錯(cuò)誤率隨迭代次數(shù)的變化曲線Fig.8 Error rate curves of broken bone segmentation (A) and broken bone splicing (B) with number of iterations

6.2.2 實(shí)驗(yàn)結(jié)果

圖9 碎骨分割與拼接結(jié)果Fig.9 Segmentation and splicing results of broken bones

本文對8例脛骨骨折臨床病例進(jìn)行分割與拼接實(shí)驗(yàn),并將該方法分割、拼接結(jié)果與骨科醫(yī)生使用的3-Matic軟件手工分割、拼接結(jié)果進(jìn)行對比,實(shí)驗(yàn)結(jié)果如圖9所示. 圖9記錄了碎骨與模板數(shù)據(jù)：算法分割、拼接結(jié)果,手工分割、拼接結(jié)果. 由圖9可見,本文方法能較好地實(shí)現(xiàn)碎骨分割與碎片拼接,接近手動分割與拼接結(jié)果.

下面對本文方法從兩方面進(jìn)行評估:

1) 本文算法的分割結(jié)果與手工分割結(jié)果進(jìn)行對比；

2) 本文算法的拼接結(jié)果與手工拼接結(jié)果進(jìn)行對比.

利用Dice系數(shù)[39]評估算法分割的精確度,實(shí)驗(yàn)結(jié)果列于表1. 表1列出了每個(gè)碎骨模型使用算法分割與手工分割兩種方法最終得到的三角面片總數(shù)量,其中算法分割三角面片為算法最終分割結(jié)果,包含分割錯(cuò)誤的三角面片. Dice系數(shù)是一種集合相似度度量函數(shù),通常用于計(jì)算兩個(gè)樣本的相似度,范圍為[0,1],Dice系數(shù)的值越接近1,說明相似度越高. 碎骨表面分割方法的Dice系數(shù)計(jì)算公式如下:

(8)

其中|Fseg∩Fgt|表示與手工分割對比本文算法正確分割三角面片個(gè)數(shù),|Fseg|表示算法分割三角面片個(gè)數(shù),|Fgt|表示手工分割三角面片個(gè)數(shù). 表1中Dice系數(shù)平均值為0.984 4. 由表1可見,本文方法分割結(jié)果與手動分割結(jié)果相比誤差較小,在可接受范圍內(nèi).

表1 碎骨分割結(jié)果

下面從平移誤差、旋轉(zhuǎn)誤差兩方面評估算法拼接的精確度,實(shí)驗(yàn)結(jié)果列于表2. 表2列出了8例碎骨所有碎片的平移誤差與旋轉(zhuǎn)誤差,并計(jì)算了平移誤差均值與旋轉(zhuǎn)誤差均值. 其中：平移誤差采用均方根誤差(RMSE)表示,定義為兩種拼接碎片對應(yīng)點(diǎn)間歐氏距離的均方根；旋轉(zhuǎn)誤差定義為算法拼接旋轉(zhuǎn)角度與手工拼接旋轉(zhuǎn)角度差值,用ICP算法計(jì)算手動拼接與算法拼接的剛體變換,先用矩陣R表示其旋轉(zhuǎn)變換,再根據(jù)旋轉(zhuǎn)矩陣R計(jì)算得到旋轉(zhuǎn)誤差α,β,θ(α,β,θ分別表示繞x,y,z軸旋轉(zhuǎn)誤差角度). 將本文方法拼接結(jié)果與手工拼接結(jié)果進(jìn)行對比,8例碎骨所有碎片間平移誤差均值為9.218 mm,旋轉(zhuǎn)平均值為1.298°,其中x,y,z軸方向上的旋轉(zhuǎn)誤差均值分別為1.515°,1.139°,1.241°. 由表2可見,本文方法對大碎片的平移、旋轉(zhuǎn)效果較好,碎片較小或碎片特征不明顯時(shí),平移、旋轉(zhuǎn)誤差較大,如實(shí)驗(yàn)1中碎片3,實(shí)驗(yàn)6中碎片3、4等. 通過與手工拼接結(jié)果比較,本文方法總體拼接誤差較小,能滿足臨床醫(yī)學(xué)碎骨拼接的要求.

表2 碎骨拼接結(jié)果

下面通過對單個(gè)骨折病例算法分割、拼接運(yùn)行時(shí)間與手工分割、拼接時(shí)間進(jìn)行對比評估算法的效率. 實(shí)驗(yàn)平臺為VSCode Linux x64,搭建在Ubuntu18.04系統(tǒng)上,算法測試環(huán)境為NVIDIA Quadro RTX 6000 GPU,系統(tǒng)內(nèi)存為16 GB. 實(shí)驗(yàn)結(jié)果列于表3. 由表3可見: 1) 本文方法對一個(gè)完整骨折病例的所有碎塊進(jìn)行分割與拼接修復(fù)時(shí)間為1～3 min,手工分割與拼接方法耗時(shí)為40～90 min,因此本文方法極大縮短了碎骨分割與拼接時(shí)間; 2) 碎骨分割、拼接耗時(shí)與碎骨碎裂個(gè)數(shù)相關(guān),骨折碎片個(gè)數(shù)越多,算法運(yùn)行時(shí)間越長. 本文方法對分割與拼接迭代進(jìn)行,結(jié)果逐步求精,碎骨分割、拼接結(jié)果與手動分割、拼接結(jié)果相近,誤差較小,耗時(shí)較短,在可接受范圍內(nèi).

表3 碎骨分割拼接時(shí)間對比

6.3 網(wǎng)格密度對比實(shí)驗(yàn)

為驗(yàn)證本文方法對不同網(wǎng)格密度的魯棒性,從高、中、低三種網(wǎng)格密度進(jìn)行驗(yàn)證,網(wǎng)格頂點(diǎn)數(shù)量分別約為20 000,10 000,5 000個(gè). 本文從Dice系數(shù)、平移誤差、旋轉(zhuǎn)誤差三方面評估不同網(wǎng)格密度下算法分割與拼接效果的魯棒性,實(shí)驗(yàn)結(jié)果列于表4.

表4 網(wǎng)格密度實(shí)驗(yàn)結(jié)果

表4列出了8例骨折數(shù)據(jù)在高、中、低3種網(wǎng)格密度下算法分割與拼接的結(jié)果,與手動分割和拼接結(jié)果進(jìn)行對比. 采用6.2.2中8例骨折數(shù)據(jù)進(jìn)行格密度對比實(shí)驗(yàn),利用在3種網(wǎng)格密度下Dice系數(shù)均值評估算法分割的魯棒性,利用平移誤差均值與旋轉(zhuǎn)誤差均值評估算法拼接的魯棒性. 由表4可見,網(wǎng)格密度越高,Dice系數(shù)越高,平移誤差與旋轉(zhuǎn)誤差越小. 圖10為實(shí)驗(yàn)過程中一例碎骨在3種不同網(wǎng)格密度下算法的分割與拼接結(jié)果.

圖10 不同網(wǎng)格密度分割與拼接結(jié)果Fig.10 Segmentation and splicing results of different grid densities

由圖10可見,網(wǎng)格密度對實(shí)驗(yàn)結(jié)果有影響,模型精度越高,算法分割與拼接效果越好. 本文算法在高、中、低3種網(wǎng)格密度下分割、拼接結(jié)果與手動分割、拼接結(jié)果相近,誤差較小,能滿足術(shù)前規(guī)劃實(shí)際要求.

6.4 碎骨拼接對比實(shí)驗(yàn)

將基于模板與數(shù)據(jù)驅(qū)動的拼接方法與ICP方法[35]以及文獻(xiàn)[40]方法進(jìn)行對比實(shí)驗(yàn). 實(shí)驗(yàn)數(shù)據(jù)采用6.2.2中8例碎骨骨折數(shù)據(jù),碎骨碎片采用手工分割方式提取. ICP方法實(shí)驗(yàn)是將碎骨碎片采用ICP算法與對側(cè)骨模板直接匹配; 文獻(xiàn)[40]方法實(shí)驗(yàn)是基于拼接思想,采用兩兩拼接方式,以兩個(gè)碎片點(diǎn)云質(zhì)心連線為掃描方向,采用距離濾波器提取兩個(gè)碎片斷裂邊界接觸區(qū)域,最終利用ICP算法匹配拼接.

下面從平移誤差與旋轉(zhuǎn)誤差兩方面評估3種拼接方法對碎骨骨折的復(fù)位效果,實(shí)驗(yàn)結(jié)果列于表5. 3種方法的平移誤差與旋轉(zhuǎn)誤差均與手工拼接方法進(jìn)行對比. 當(dāng)碎骨碎片之間相對旋轉(zhuǎn)角度較大或者碎片間距較大時(shí),ICP方法與文獻(xiàn)[40]方法都對初始位置較敏感,需進(jìn)行手工粗匹配,圖11為兩例骨折數(shù)據(jù)在3種拼接方法下的實(shí)驗(yàn)結(jié)果. ICP方法與文獻(xiàn)[40]方法對于多個(gè)碎片拼接的平移誤差與旋轉(zhuǎn)誤差均較大,對小碎片位置匹配較差; 本文方法利用碎骨與模板間的特征密集對應(yīng)關(guān)系進(jìn)行粗匹配,再根據(jù)ICP算法進(jìn)行精確匹配,碎片匹配精確度較高,拼接誤差較小,可滿足骨折術(shù)前規(guī)劃的要求.

表5 不同拼接方法實(shí)驗(yàn)結(jié)果對比

圖11 3種方法拼接結(jié)果對比Fig.11 Comparison of splicing results of 3 methods

綜上所述,針對骨折發(fā)生隨機(jī)、相互接觸碎骨的邊界分割困難、小碎骨丟失導(dǎo)致碎骨拼接較難的問題,本文提出了一種基于模板引導(dǎo)與數(shù)據(jù)驅(qū)動的三維碎骨修復(fù)方法. 該方法的核心思想是基于數(shù)據(jù)驅(qū)動初步提取碎骨表面,利用特征將碎骨與模板進(jìn)行匹配,通過分割-拼接迭代求精,實(shí)現(xiàn)碎骨提取與碎骨拼接. 用8例脛骨骨折臨床數(shù)據(jù)進(jìn)行驗(yàn)證,所有實(shí)驗(yàn)結(jié)果均符合骨折修復(fù)術(shù)前規(guī)劃實(shí)際要求,與手工分割與拼接結(jié)果相比,誤差較小.

與目前已有方法相比,本文方法優(yōu)點(diǎn)如下:

1) 實(shí)現(xiàn)了碎骨表面模型的提取,并完成了碎骨拼接修復(fù)；

2) 基于模板匹配更好地實(shí)現(xiàn)了碎骨邊界分割,解決了因碎骨邊界接觸、邊界特征不明顯導(dǎo)致邊界分割困難的問題；

3) 基于模板匹配實(shí)現(xiàn)碎骨修復(fù)可更好地解決因細(xì)小碎片丟失導(dǎo)致某些碎片邊界吻合程度不高、碎片拼接困難的問題.

但本文方法還存在以下不足: 由于骨折發(fā)生隨機(jī),存在碎骨碎片較小,碎片特征不明顯,導(dǎo)致小碎片特征匹配精度不高,匹配誤差較大,所以小碎片的分割與拼接誤差相對較大.