羅月童，徐冰清，薛 浩，石放放，張 偉，丁 勝

(1. 合肥工業(yè)大學計算機與信息學院，安徽 合肥 230009；2. 合肥工業(yè)大學信息工程系，安徽 宣城 242099；3. 安徽醫(yī)科大學第四附屬醫(yī)院，安徽 合肥 230022)

基于對稱性的頜骨復位方法研究及應(yīng)用

羅月童1,2，徐冰清1，薛 浩3，石放放1，張 偉1，丁 勝1

(1. 合肥工業(yè)大學計算機與信息學院，安徽 合肥 230009；2. 合肥工業(yè)大學信息工程系，安徽 宣城 242099；3. 安徽醫(yī)科大學第四附屬醫(yī)院，安徽 合肥 230022)

針對嚴重錯位骨折的頜骨復位口腔外科手術(shù)，提出一種幫助醫(yī)生進行術(shù)前規(guī)劃和準備的破損頜骨復位方法。首先使用移動立方體(MC)算法從患者CT圖像中提取頜骨的面片模型；然后結(jié)合圖割算法、形狀直徑函數(shù)(SDF)描述子和混合高斯模型的交互式分割方法，允許醫(yī)生通過簡單的勾畫快速分割出破損骨塊；最后基于頜骨的對稱性特性自動完成破損骨塊復位。通過破損骨塊復位完成頜骨模型重構(gòu)，能完整保留對醫(yī)生非常重要的骨折線信息。該方法已應(yīng)用于自主開發(fā)的計算機輔助破損頜骨復位系統(tǒng)，并使用安徽醫(yī)科大學附屬醫(yī)院的患者數(shù)據(jù)進行了測試，結(jié)果表明了該方法的有效性。

計算機輔助手術(shù)；頜骨復位；網(wǎng)格分割；對稱性檢測

頜骨(如圖1所示)因具有形態(tài)不規(guī)則、骨縫多、骨質(zhì)薄弱等特點，易在運動、交通事故等過程中受到撞擊而出現(xiàn)嚴重移位骨折，因此頜骨的復位與固定是非常普遍的口腔科手術(shù)[1]。CT圖像已被廣泛用于頜骨外科手術(shù)的術(shù)前檢查和術(shù)后效果評估，近年來計算機輔助手術(shù)(computer aided surgery，CAS)領(lǐng)域探索基于 CT數(shù)據(jù)進一步幫助醫(yī)生進行手術(shù)規(guī)劃、術(shù)前準備等的方法，以進一步減少頜骨手術(shù)時間、提高手術(shù)質(zhì)量[2]。

圖1 頜骨示意圖

Patel等[3]較早地研究了基于CT圖像的頜骨手術(shù)規(guī)劃和模擬系統(tǒng)的關(guān)鍵技術(shù)；Cevidanes等[4]研究了下頜骨手術(shù)中的變形計算中相關(guān)的圖像處理技術(shù)；Enciso等[5]和Mollemans等[6]都針對咬合運動模擬中的軟組織建模問題開展了研究；Chowdhury等[7-8]基于患者 CT圖像為患者重建完整的正常頜骨模型，輔助醫(yī)生進行術(shù)前規(guī)劃。目前，因受時間、場地等因素的影響，多數(shù)臨床醫(yī)生沒能使用手術(shù)模擬系統(tǒng)，但是完整的患者頜骨模型能方便地被醫(yī)生采用，尤其隨著3D打印技術(shù)的普及，完整的頜骨模型可幫助醫(yī)生進行術(shù)前規(guī)劃、術(shù)前準備及與患者進行溝通[9]。本文針對快速、方便重構(gòu)頜骨模型的關(guān)鍵技術(shù)開展了研究。

因為患者兩側(cè)頜骨很少同時受傷且正常人的頜骨具有對稱性，因此利用鏡像原理重構(gòu)完整頜骨模型是常用的方法[9]，但重構(gòu)所得模型缺少骨折面信息，而骨折面對醫(yī)生選擇定位位置、預彎鈦板等有重要參考價值。文獻[7-8]首先對患者CT圖像進行精確分割，然后依據(jù)分割圖像重構(gòu)骨折面，最后通過使用迭代最近點(iterative closest point，ICP)算法匹配骨折面，將錯位骨塊恢復至正常位置，從而實現(xiàn)完整頜骨模型的重構(gòu)。其方法的效果依賴2點：①CT圖像的準確分割，保證準確構(gòu)建骨折面；②骨折面的準確匹配。但CT圖像的自動準確分割還是一個沒有完全解決的問題，而且因頜骨較薄，骨折面較窄，對很窄的面進行匹配可能帶來誤差，因此文獻[7-8]的方法難以保證精度，或需要較多的后續(xù)處理才能滿足要求。本文方法基于頜骨的對稱性恢復錯位骨塊的位置，對稱信息屬于頜骨的宏觀信息，具有較好的魯棒性。

1 方法框架

本文方法的整體流程如圖 2所示，即：①讀入患者的CT圖像數(shù)據(jù)。因為對于頜骨修復只需要關(guān)注其外觀信息，無需處理骨骼內(nèi)部數(shù)據(jù)，因此本文使用移動立方體算法(marching cube，MC)將CT圖像數(shù)據(jù)轉(zhuǎn)換成頜骨的面片模型，并使用邊折疊算法[10]對 MC算法所提取模型進行優(yōu)化，生成最終的頜骨面片模型。②對頜骨面片模型進行破損骨塊分割，獲取一組骨塊。③基于對稱性原理將錯位骨塊恢復至正常位置，從而構(gòu)建完整的頜骨模型。④輸出完整的頜骨模型。整個過程是一個交互過程，用戶如果觀察到分割結(jié)果或重構(gòu)結(jié)果不符合要求，可以對分割過程或復位過程進行手工干預。

圖2 本文方法整體流程

2 基于筆觸的半自動骨塊分割

2.1 方法選擇

由圖2可知，骨塊分割的本質(zhì)是網(wǎng)格分割。網(wǎng)格分割因為一直是圖形學領(lǐng)域的重要問題而備受關(guān)注[11]。早期的網(wǎng)格分割算法多為手工分割，但手工分割費時費力，難以被非專業(yè)用戶所接受。雖然人們對自動網(wǎng)格分割算法開展了大量研究，但對復雜網(wǎng)格模型進行完全自動分割仍是一個挑戰(zhàn)[11]，因此交互式半自動網(wǎng)格分割算法最近頗受關(guān)注[12]，基于筆觸的半自動分割方法就是其中最重要的方法之一[13]，因此本文采用基于筆觸的半自動分割方法實現(xiàn)骨塊的快速分割。

基于筆觸的分割方法由用戶通過繪制指定部分前景網(wǎng)格和背景網(wǎng)格，然后基于某種算法實現(xiàn)整個網(wǎng)格的自動分割，其中有 2個關(guān)鍵點：①如何刻畫網(wǎng)格的局部特征，經(jīng)常被采用的局部信息包括基于曲率圖的測地距離[14]、Isophotic度量距離[15]、Lai等[16]提出的特征敏感度量、diffusion distance[17]、參數(shù)曲面的各向異性測地距離[18]、形狀直徑函數(shù)(shape diameter function，SDF)[19]等。②如何將用戶的少量輸入擴展到整個模型，采用算法包括區(qū)域增長[20-21]、圖割(graph-cut)[22]、random walks[23-24]、hierarchical aggregation[25]等。因為破損的骨塊一定是在連接處發(fā)生了突變，所以本文選擇如下：

(1) 選用 SDF作為網(wǎng)格的局部特征。SDF值能將模型的體積信息反映到表面模型上，破損的骨塊在連接處的形狀突變必然引起網(wǎng)格SDF值的變化，此變化將成為模型分割的依據(jù)。

(2) 選擇graph-cut算法。graph-cut通過最大流/最小割算法思想解決圖的平滑切割問題，可以綜合考慮網(wǎng)格自身的局部信息及相鄰網(wǎng)格之間的連接關(guān)系。分割破損骨塊的主要依據(jù)是其連接處的突變，因此將相鄰網(wǎng)格的連接關(guān)系納入考慮，符合骨塊分割的特點。

2.2 算法步驟

文獻[13]所提出的分割方法和本文的需求比較吻合，因此本文基于文獻[13]的方法實現(xiàn)骨塊分割算法，具體過程如下：

步驟1. 計算SDF值。計算頜骨表面模型S中每個三角面片 fi(fi∈S)的SDF值，記為 SDF（fi）。

步驟2. 筆觸交互。用戶使用鼠標交互地在需要被分割的破損骨塊上繪制若干筆觸，假設(shè)和筆觸相交的所有三角面片屬于需要被分割骨塊，并記為F。

步驟 3. 建立前景高斯混合模型。對所有fj∈F的 SDF（fj）使用極大似然估計法建立高斯混合模型，記為前景高斯混合模型Gf。Gf的概率分布密度函數(shù)為：

其中，

表示第j個分量的概率分布函數(shù)；本文方法前景高斯混合模型有2個分量，即K=2。

步驟4. 建立背景高斯混合模型。和步驟3類似，對所有 fj∈（S -F）的 SDF（fj）使用極大似然估計法建立高斯混合模型，記為背景高斯混合模型Gb。Gb的概率分布密度函數(shù)為：

其中，背景高斯混合模型有4個分量，即K=4。

步驟5. 使用graph-cut算法分割。對表面模型S進行“前景-背景”分割，分割所得前景即為所求破損骨塊。graph-cut算法優(yōu)化式(3)所示能量函數(shù)最小：

其中，lf=1表示三角面片f被標記為前景，lf=0表示三角面片f被標記為背景。

其中，θ（fp,fq）表示相鄰面片fp和fp之間的夾角。

步驟6. 觀察分割結(jié)果，以圖3為例。圖3(a)為未分割模型，紅色虛線內(nèi)為骨折處。如果分割出的結(jié)果不完整或有錯誤，如圖3(b)所示，則返回步驟2繼續(xù)優(yōu)化筆觸；如果分割出的結(jié)果已經(jīng)完全包含破損骨塊，如圖3(c)所示，則保留此分割結(jié)果，在此基礎(chǔ)上進行分割，返回步驟2進行筆觸交互；如果分割正確，如圖3(d)所示，則轉(zhuǎn)向步驟7。

步驟7. 保存分割所得破損骨塊，如果還有破損骨塊需要分割，則轉(zhuǎn)向步驟2，否則結(jié)束。

圖3 分割情況示意圖

3 基于對稱性的骨塊復位

因為頜骨具有對稱性，且患者很少出現(xiàn)兩側(cè)頜骨同時受損的情況，因此可以通過找到破損骨塊的對稱骨塊，然后基于對稱骨塊將破損骨塊復位至正常位置，這涉及 2個問題：①如何找到破損骨塊對應(yīng)的正常骨塊，即對稱性檢測問題；②如何基于正常骨塊將破損骨塊恢復至正常位置。

文獻[26]詳細總結(jié)了對稱性檢測，文獻[27]則通過處理采樣點對之間的對稱矩陣，檢測出模型中任意全局和局部對稱性，同時計算出對稱部分之間的變換矩陣，不僅檢測破損骨塊對應(yīng)的正常骨塊，而且可以利用變換矩陣直接實現(xiàn)破損骨塊的復位(圖 4)。本文方法基于文獻[27]實現(xiàn)破損骨塊復位，具體過程如下：

步驟1. 預處理。

(1) 采樣：從頜骨表面模型S中隨機抽取N個頂點作為樣本點，用SP表示樣本點集。綜合考慮計算速度和精度，N通常取表面模型S頂點總數(shù)的1/10。

圖4 利用對稱性復位的具體過程

(2) 特征計算：

①構(gòu)造特征矩陣：為每個pi∈SP按式(6)構(gòu)造4×4的特征矩陣 E (pi)：

其中，B表示以樣本點pi為中心，以一個給定值為半徑的球，表示該球的表面積，e為 B中網(wǎng)格的邊，為e方向上的單位向量，為e∩B的長度，β(e)是以e為公共邊的兩個三角形法向的夾角的補角。

②構(gòu)造樣本點的局部坐標系：若特征矩陣E (pi)是非滿秩矩陣，則認為該樣本點為無效樣本，進行剔除，否則求解該矩陣的特征值和特征向量，并將最小特征值對應(yīng)的特征向量作為局部坐標系的 Z軸，將剩余兩個特征值對應(yīng)的特征向量分別作為X軸和Y軸。另外，將最小特征值和最大特征值分別作為采樣點pi的最小曲率rmin和最大曲率rmax。若，則該采樣點無效，其中0＜ε≤1，ε默認值為0.75。

步驟2. 對整個頜骨模型，按下列步驟獲得主對稱面信息F。

(1) 從樣本點集SP中挑出所有屬于主骨塊M的樣本點，記為SPM，。

(3) 使用多維尺度(multi dimensional scaling，MDS)算法將所有的 Q(pi,pj)降維到二維空間，并繪制散點圖；

(4) 用戶在散點圖上交互勾畫出點聚集密度最高的區(qū)域，用該區(qū)域中所有的點的平均值作為整個頜骨模型的主對稱信息(A、B、C、D)，用該對稱信息構(gòu)造矩陣：

步驟3. 對破損骨塊pi，按下列步驟獲得將其復位到正確位置的變換矩陣Ti。

(1)從樣本點集SP中挑出所有屬于破損骨塊pi的樣本點，記為SPPi，。

(2) 從樣本點集SP中挑出所有屬于主骨塊M的樣本點，記為SPM，。

(4) 使用 MDS算法將所有的 M（pi,pj）降維到二維空間，并繪制散點圖。

(5) 用戶在散點圖上交互勾畫出點密度最高的區(qū)域，用該區(qū)域中所有的點的平均值作為破損骨塊pi的變換矩陣Ti。

4 應(yīng) 用

本文基于 Voreen實現(xiàn)頜骨復位算法，Voreen是一個開放源碼的跨平臺可視化編程庫，由明斯特大學計算機系的可視化與計算機圖形學研究中心使用 C++語言開發(fā)。Voreen提供現(xiàn)有的體繪制算法，也為研究者提供開發(fā)新算法的接口，使研究者能將自己的算法集成到Voreen框架中[28]。本系統(tǒng)的開發(fā)環(huán)境為Visual Studio 2010，在配置為Intel2.93 GHz的CPU以及NVidia GeForce GTS450顯卡的個人電腦上能流暢運行。

本文使用來自安徽醫(yī)科大學第一附屬醫(yī)院的真實數(shù)據(jù)進行測試。圖5(a)為原始模型，其下頜骨中部錯位，圖5(b)是分割結(jié)果；圖5(c)是復位后的頜骨模型。

圖5 實驗效果圖

5 結(jié) 論

計算機輔助缺損頜骨修復技術(shù)能完成復雜骨折的個體化復位和解剖重建，對面中部骨折的治療具有非常實用的臨床價值，同時其所涉及的方面也是計算機視覺中的一個突出難點。在計算機視覺方面，本文實現(xiàn)了一種基于筆觸的網(wǎng)格分割方法，對基于對稱性的頜骨復位技術(shù)進行了進一步研究；CAS方面，本文實現(xiàn)了一個計算機輔助破損頜骨復位系統(tǒng)，對頜面樣本數(shù)據(jù)處理效果良好，可以輔助醫(yī)生提前做好術(shù)前準備和溝通工作。由于頜骨骨折的情況因個體情況不同差異較大，本文系統(tǒng)仍需要在實際應(yīng)用中不斷完善與優(yōu)化，提供更準確、全面的輔助信息。

[1] King R E, Scianna J M, Petruzzelli G J. Mandible fracture patterns: a suburban trauma center experience [J]. American Journal of Otolaryngology, 2004, 25(5): 301-307.

[2] Chowdhury A S, Bhandarkar S M. Computer vision-guided virtual craniofacial surgery: a graph-theoretic and statistical perspective [M]. London: Springer-Verlag London Limited, 2011: 1-2.

[3] Patel V V, Vannier M W, Marsh J L, et al. Assessing craniofacial surgical simulation [J]. Computer Graphics & Applications IEEE, 1996, 16(1): 46-54.

[4] Cevidanes L H S, Styner M A, Phillips C, et al. 3D morphometric changes 1 year after jaw surgery [C]// Biomedical Imaging: from Nano to Macro, 2007. ISBI 2007. 4th IEEE International Symposium on. New York: IEEE Press, 2007: 1332-1335.

[5] Enciso R, Memon A, Neumann U, et al. The virtual craniofacial patient project: 3D jaw modelling and animation [J]. Studies in Health Technology and Informatics, 2003, 94: 65-71.

[6] Mollemans W, Schutyser F, Nadjm i N, et al. Very fast soft tissue predictions with mass tensor model for maxillofacial surgery planning systems [J]. International Congress, 2005, 1281: 491-496.

[7] Chowdhury A S, Bhandarkar S M. Virtual craniofacial reconstruction using computer vision, graph theory and geometric constraints [J]. Pattern Recognition Letters, 2009, 30(10): 931-938.

[8] Bhandarkar S M, Chowdhury A S, Tang Y. Computer vision guided virtual craniofacial reconstruction [J]. Computerized Medical Imaging & Graphics, 2007, 31(6): 418-427.

[9] Liu H Q, Hu Q X, Li L, et al. A study of the method of reconstructing the bionic scaffold for repairing defective bone based on tissue engineering [C]//Proceedings of PROLAMAT 2006, IFIP TC5 International Conference, June 15-17, 2006, Shanghai, China. Boston: Springer US, 2006: 650-657.

[10] Hoppe H. Progressive meshes [C]//SIGGRAPH’96 Proceeding of the 23rdAnnual Conference on Computer Graphics and Interactive Techniques. New York: ACM Press, 1996: 99-108.

[11] Sham ir A. A survey on mesh segmentation techniques [J]. Computer Graphics Forum, 2008, 27(6): 1539-1556.

[12] 秦緒佳, 歐宗瑛, 紀鳳欣, 等. 醫(yī)學圖像的交互分割及三維表面重建[J]. 工程圖學學報, 2001, 22(2): 94-101.

[13] Fan L B, Li L G, Liu K. Paint mesh cutting [J]. Computer Graphics Forum, 2011, 30(2): 603-612.

[14] Gatzke T, Grimm C, Garland M, et al. Curvature maps for local shape comparison [C]//Shape Modeling and Applications, 2005 International Conference. Los A lam itos: IEEE Computer Society Press, 2005: 244-253.

[15] Pottmann H, Steiner T, Hofer M, et al. The isophotic metric and its application to feature sensitive morphology on surfaces [J]. Lecture Notes in Computer Science, 2004, 3024: 560-572.

[16] Lai Y K, Zhou Q Y, Hu S M, et al. Robust feature classification and editing [J]. IEEE Transactions on Visualization & Computer Graphics, 2007, 13(1): 34-45.

[17] de Goes F, Goldenstein S, Velho L. A hierarchical segmentation of articulated bodies [J]. Symposium on Geometry Processing, 2008, 27(5): 1349-1356.

[18] Seong J K, Jeong W K, Cohen E. Anisotropic geodesic distance computation for parametric surfaces [C]//Shape Modeling and Applications, 2008. SM I 2008. IEEE International Conference on. New York: IEEE Press, 2008: 179-186.

[19] Shapira L, Sham ir A, Cohen-Or D. Consistent mesh partitioning and skeletonisation using the shape diameter function [J]. Visual Computer, 2008, 24(4): 249-259.

[20] Ji Z P, Liu L G, Chen Z G, et al. Easy mesh cutting[J]. Computer Graphics Forum. 2006, 25(3): 283-291.

[21] Wu H Y, Pan C H, Pan J, et al. A sketch-based interactive framework for real-time mesh segmentation [EB/OL]. [2015-10-15]. http://nlpr-web.ia.ac.cn/2007papers/gjhy/ gh108.pdf.

[22] Brown S, Morse B, Barrett W. Interactive part selection for mesh and point models using hierarchical Graph-cut partitioning [C]//Proceedings of Graphics Interface. New York: ACM Press, 2009: 23-30.

[23] Lai Y K, Hu S M, Martin R R, et al. Fast mesh segmentation using random walks [C]//Proceedings of the 2008 ACM Symposium on Solid and Physical Modeling. New York: ACM Press, 2008: 183-191.

[24] Zhang J Y, Wu C L, Cai J F, et al. Mesh snapping: robust interactive mesh cutting using fast geodesic curvature flow [J]. Computer Graphics Forum, 2010, 29(2): 517-526.

[25] Xiao C X, Fu H B, Tai C L. Hierarchical aggregation for efficient shape extraction [J]. The Visual Computer, 2009, 25(3): 267-278.

[26] 史 川, 葉志偉, 徐鵬捷. 對稱性檢測方法的研究分析[J]. 現(xiàn)代計算機: 專業(yè)版, 2010, 9(9): 67-69.

[27] Mitra N J, Guibas L J, Pauly M. Partial and approximate symmetry detection for 3D geometry [J]. ACM Transactions on Graphics, 2006, 25(3): 560-568.

[28] Meyer-Spradow J, Ropinski T, Mensmann J, et al. Voreen: a rapid-prototyping environment for ray-casting-based volume visualizations [J]. IEEE Computer Graphics and Applications, 2009, 29(6): 6-13.

Symmetry Based Broken Jaw Reconstruction Method and its App lication

Luo Yuetong1,2, Xu Bingqing1, Xue Hao3, Shi Fangfang1, Zhang Wei1, Ding Sheng1

(1. School of Computer and Information, Hefei University of Technology, Hefei Anhui 230009, China; 2. Department of Information Engineering, Hefei University of Technology, Xuancheng Anhui 242099, China; 3. The Fourth A ffiliated Hospital of Anhui Medical University, Hefei Anhui 230022, China)

For oral surgery of dislocated fracture jaw reconstruction, a broken jaw reconstruction method is presented to help doctors for preoperative planning and preparation. Firstly, marching cube algorithm is used to extract facet model from patient’s CT image data, then a sketch-based interactive segmentation method, which combines graph cut algorithm and SDF descriptor and Gaussian m ixture models, is presented to allow doctors to segment broken jaw rapidly by simple strokes; at last reconstruction is realized automatically based on the craniofacial symmetry. This method accomplishes reconstruction by reposition of broken jaw so that it can fully keep important fracture information for doctors. The method has been applied in self-developed jaw reconstruction system; patient data from The First A ffiliated Hospital of Anhui Medical University has been used for testing, the result demonstrates the effectiveness of this approach.

computer aided surgery; jaw reconstruction; mesh segmentation; symmetry detection

TP 391

10.11996/JG.j.2095-302X.2016030316

A

2095-302X(2016)03-0316-06

2015-11-06；定稿日期：2015-11-25

國家自然科學基金項目(11305205，61370167，61305093)；中石油-中國科學院重大戰(zhàn)略合作項目(2015A-4812)

羅月童(1978–)，男，安徽青陽人，教授，博士。主要研究方向為科學計算可視化、計算機圖形學和計算機輔助設(shè)計。E-mail：ytluo@hfut.edu.cn