馬 天,楊 秦,李占利
(西安科技大學(xué)計算機科學(xué)與技術(shù)學(xué)院,陜西 西安 710054)
近年來,隨著計算機輔助設(shè)計仿真、三維重建技術(shù)與快速成型技術(shù)的快速發(fā)展,無托槽隱形矯治引起了廣大正畸醫(yī)生以及患者的格外關(guān)注,其具有美觀、舒適、安全、便捷等優(yōu)點[1]。虛擬牙齒正畸系統(tǒng)實現(xiàn)了在三維可視化狀態(tài)下牙齒移動的方式和步驟的設(shè)計,主要是為制作隱形矯治器提供數(shù)據(jù)源[2]。該系統(tǒng)首先通過先進的掃描設(shè)備獲取患者牙齒的牙頜數(shù)字模型[3],其次在虛擬牙齒正畸系統(tǒng)中進行單顆牙齒分割修補[4]、牙齒移動和矯治方案制定、牙齦組織變形等[5]環(huán)節(jié),最終輸出隱形矯治器的母版。其中牙齒移動和矯治方案的確定是非常重要的技術(shù)環(huán)節(jié),通過計算機算法實現(xiàn)正畸中間路徑的獲取,可以減少對口腔醫(yī)師經(jīng)驗的依賴,提高正畸的精準度,因此牙齒矯治自動化是口腔醫(yī)學(xué)發(fā)展的趨勢[6]。
牙齒正畸是三維環(huán)境下有障礙物的多目標路徑規(guī)劃問題[7],傳統(tǒng)的路徑規(guī)劃優(yōu)化方法在應(yīng)對復(fù)雜環(huán)境時會存在一定的缺陷[8]。在處理復(fù)雜環(huán)境信息路徑規(guī)劃時,來自于自然界的啟示往往能起到很好的作用,常用的仿生學(xué)算法有:蟻群算法、神經(jīng)網(wǎng)絡(luò)算法、粒子群算法和遺傳算法等,其中粒子群采用實數(shù)編碼,更方便問題的映射。
國內(nèi)外有關(guān)自動牙齒正畸路徑規(guī)劃的研究成果較少,MOTOHASHI 和KURODA[9]的方法首先計算每顆牙齒的特征線段,通過牙弓線牽引的方式實現(xiàn)自動排牙技術(shù),該方法不適宜牙列擁擠碰撞的情況,因此在實際正畸中存在諸多問題;LIN[10]提出半自動排牙技術(shù),仍需提取單顆牙齒特征線段,利用二次函數(shù)建立擬合牙弓線,未規(guī)避牙齒間可能產(chǎn)生的碰撞和牽引移動過程中牙齒的生理特征,采用該正畸方案可能會帶來負面效果。從1998 年開始,國外開始出現(xiàn)隱形矯治器虛擬牙齒正畸系統(tǒng),ALIGN 公司[11]自主研發(fā)的Invisalign 在正畸臨床上取得了巨大成功;1999 年,GeoDigm 公司[12]開發(fā)出一款可以操作牙齒分割、移動及旋轉(zhuǎn)的虛擬牙齒正畸軟件Emodel;2003 年,美國Cadent 公司[13]實現(xiàn)了牙齒正畸參數(shù)分析和牙齒正畸自動排牙算法。國外關(guān)于虛擬牙齒正畸系統(tǒng)的文獻大多為公司相關(guān)產(chǎn)品的宣傳文案,而有關(guān)虛擬牙齒正畸系統(tǒng)技術(shù)層面的文獻很少。
近年來,國內(nèi)在牙齒正畸方面也取得了一些成果,2010 年LI 等[14]采取遺傳算法求解牙齒正畸路徑規(guī)劃,認為遺傳算法支持多個物體同時移動符合牙齒正畸的需求,但遺傳算法存在多參數(shù)的自適應(yīng)性問題且效率較低;2016 年山東大學(xué)張筱[15]提出單顆牙齒的排牙方法和牙齒正畸過程中的碰撞檢測方法,并設(shè)計牙齒正畸路徑規(guī)劃的整體思路流程;2018 年付敬鼎[16]提出基于改進快速搜索隨機樹(rapidly-exploring random trees,RRT)算法進行路徑規(guī)劃,由于算法本身的隨機性,生成的路徑比較曲折,甚至出現(xiàn)繞遠路。2020 年徐曉強等[17]采用粒子群的算法,以一個粒子代表整口牙齒進行矯治,存在維度過高導(dǎo)致運算效率較低的問題。
根據(jù)臨床經(jīng)驗,牙齒正畸需要考慮其生理組織結(jié)構(gòu),在以往的研究中,未對不同種類的牙齒進行區(qū)分,導(dǎo)致實際和理想正畸過程存在一定的差異。本文改進多粒子群牙齒正畸路徑規(guī)劃的方法考慮了多顆牙齒正畸的高維度和碰撞問題,將每顆牙齒采用一個粒子群的方式進行矯治來降低維度以提高算法的效率;并考慮到牙齒的生理組織結(jié)構(gòu),改進不同粒子群的慣性權(quán)值和不同正畸階段位置更新的上下限,也達到了減少碰撞和優(yōu)化正畸效果的目的。
無托槽隱形矯治需要解決牙齒理想位姿確定、牙齒正畸路徑規(guī)劃、牙齦組織變形等諸多技術(shù)難題,涉及計算機圖形學(xué)領(lǐng)域的諸多分支,是一個多學(xué)科相結(jié)合的復(fù)雜系統(tǒng)。本課題主要研究的是在牙齒初始位置和目標位置均已確定的情況下,獲取牙齒正畸中間階段的位置。因此圍繞牙齒矯治過程,將解決實際牙齒三維環(huán)境下有障礙的路徑規(guī)劃為本文的研究重點。
牙齒正畸路徑規(guī)劃是獲取牙齒從初始位置到理想位置的過程,因此確定理想位置的排牙技術(shù)是無托槽隱形矯治技術(shù)的首要前提。牙齒排列是指在牙齒初始位置已知的情況下,參照排牙原則及理想的牙弓曲線[18],經(jīng)過排牙處理得到牙齒的最終位置。
牙齒的理想位置是矯治的重要參考,正確地確定患者牙弓形態(tài)極為重要,其是建立穩(wěn)定、緊密咬合關(guān)系的重要基礎(chǔ)。許多專家研究了人類牙弓線數(shù)學(xué)模型,觀點不一,主要用到的數(shù)學(xué)模型有橢圓線函數(shù)、拋物線函數(shù)、多項式函數(shù)、Beta方程等[19]。本文選用Beta方程作為擬合函數(shù)模型,由上下頜牙齒特征點集在橫斷面投影點擬合曲線,計算各投影點到曲線的最短距離得到平移量Δx,Δy;平移量Δz,則是將牙齒特征點集投影至矢狀面,并通過擬合計算得到,其反映的Spee 曲線[20]的走勢;最終通過計算的平移量來確定理想位置,圖1 為牙齒排列的曲線擬合過程。
圖1 牙齒排列技術(shù)((a)下頜牙齒橫斷面;(b)下頜牙齒橫斷面投影點擬合曲線;(c)下頜牙齒矢狀面;(d)下頜牙齒矢狀面投影點擬合曲線) Fig.1 Tooth arrangement technique ((a) Cross-section of the mandibular teeth;(b) Projection point fitting curve of cross section of mandibular teeth;(c) Sagittal plane of mandibular teeth;(d) Fitting curve of projection points on sagittal plane of mandibular teeth)
牙齒的旋轉(zhuǎn)值在口腔醫(yī)學(xué)中被稱為轉(zhuǎn)矩角、軸傾角、旋轉(zhuǎn)角[21],通過牙齒在空間中任意姿態(tài)的變換可由歐拉角構(gòu)造旋轉(zhuǎn)變換矩陣求出3 個旋轉(zhuǎn)量,圖2 為牙齒旋轉(zhuǎn)示意圖。
圖2 牙齒旋轉(zhuǎn)分量示意圖((a)歐拉角坐標;(b)牙齒旋轉(zhuǎn)坐標) Fig.2 Schematic diagram of tooth rotation component ((a) Euler angular coordinates;(b) Spatial coordinates of tooth rotation)
牙齒碰撞是指多顆牙齒在正畸路徑上發(fā)生交叉時可能出現(xiàn)的情形,牙齒正畸治療的一大難點是如何避免牙齒間的碰撞。由于牙齒本身的不規(guī)則性,且其使用的都是大數(shù)據(jù)量的三角網(wǎng)格,碰撞檢測需要消耗大量時間,因此需要對正畸方案中的碰撞檢測進行調(diào)整和簡化。
基于包圍盒的碰撞檢測算法是一類重要的檢測算法[22],其原理是在三維圖像場景中利用簡單的幾何體即包圍盒包圍碰撞檢測的幾何對象,然后將包圍盒進行碰撞檢測的結(jié)果用以替代幾何對象的碰撞檢測結(jié)果。包圍盒的形狀越簡單,碰撞檢測算法的效率越高,若包圍盒形狀越接近于研究幾何對象,其越接近真實幾何對象的碰撞檢測結(jié)果,所以應(yīng)盡量選擇形狀簡單且接近幾何對象的包圍盒。常見的包圍盒算法有軸對齊包圍盒(axis-aligned bounding box,AABB)、包圍球、方向包圍盒(oriented bounding box,OBB)以及固定方向凸包(fixed directions hulls,F(xiàn)DH)[23]。牙齒正畸由于空間較為狹窄,需選用盡可能緊密的包圍盒十分必要,且OBB 相交檢測時只需要檢測15 條分離軸,算法復(fù)雜度低,因此本文采用OBB 進行碰撞檢測,圖3為采用MATLAB 創(chuàng)建的OBB。
圖3 OBB 創(chuàng)建圖 Fig.3 OBB creation diagram
對牙齒進行矯治時,需要考慮如生理組織結(jié)構(gòu)重建、牙齒能承受的正畸力范圍、個體因素等,由于這些條件的限制,牙齒在正畸時需要分階段進行,并且每個正畸階段的移動量有限。在理想的正畸階段中,各顆牙齒移動時應(yīng)該感受不到相鄰牙齒的存在,各自朝著預(yù)定的路徑移動,不會發(fā)生碰撞。
設(shè)牙齒為Ti(i=1,2,…,n),患者牙齒個數(shù)為n,矯治階段為m,牙齒正畸是牙齒從初始姿態(tài)Pi0到目標姿態(tài)Pim的碰撞路徑,牙齒在第m階段的旋轉(zhuǎn)中心為Ci(m)=(Cxi(m),Cyi(m),Czi(m)),姿態(tài)角為δim(m)=(αi(m),βi(m),γi(m)),牙齒Ti的第m個正畸階段位姿可以表示為Pim=(Ci(m),δi(m)),則牙齒Ti第m個正畸階段的路徑規(guī)劃長度為
為了簡化牙齒旋轉(zhuǎn)的約束條件,則牙齒Ti第m個正畸階段旋轉(zhuǎn)角度為
根據(jù)臨床正畸經(jīng)驗,牙齒在正畸的一個階段中產(chǎn)生的平移或旋轉(zhuǎn)的角度不能超過一定范圍,單階段牙齒平移的值Sim不能超過0.5,δim不能超過3°。
牙齒在正畸過程中,若正畸方案不當,可能使牙齒與其相鄰牙齒發(fā)生碰撞干涉,影響正畸效果。因此考慮碰撞約束,能夠確保牙齒沿正畸路徑移動時不與相鄰牙齒發(fā)生碰撞干涉,設(shè)Kim為判斷牙齒Ti在第m階段是否發(fā)生碰撞
為了讓相鄰牙齒發(fā)生碰撞時,碰撞罰函數(shù)Fim能夠有明顯的變化,則給定一個較大常數(shù)項L,發(fā)生碰撞時罰函數(shù)為
在牙齒正畸時,牙齒之間間隙很小,相應(yīng)的活動空間也很小,因此牙齒正畸路徑規(guī)劃問題的優(yōu)化指標有:所有牙齒移動路徑最短、旋轉(zhuǎn)角度最小等;且需滿足式(4)的碰撞約束。由式(1)和(2)可得,n顆牙齒m個階段牙齒的最短移動路徑和最小旋轉(zhuǎn)角度分別為
粒子群算法初始化為一群隨機的粒子(隨機解),根據(jù)迭代找到最優(yōu)解。所有的粒子都有一個由被優(yōu)化的函數(shù)決定的適應(yīng)值,每個粒子還有一個速度決定其飛出的方向和距離,然后粒子就追隨當前的最優(yōu)粒子在解空間中搜索。粒子具有速度和位置2 個屬性,速度代表移動的快慢,位置代表移動的方向。每個粒子單獨搜尋的最優(yōu)解叫做個體極值,粒子群中最優(yōu)的個體極值作為當前全局最優(yōu)解,不斷迭代更新速度和位置,最終得到滿足終止條件的最優(yōu)解。
假設(shè)在一個D 維搜索空間中,有Z個粒子組成一個群落,其中第j個粒子表示為一個D 維的向量,則粒子j當前位置為Xj=(xj1,xj2,…,xjD),當前速度為Vj=(vj1,vj2,…,vjD),該粒子j經(jīng)歷的最優(yōu)位置為Pj=(pj1,pj2,…,pjD),該粒子群的最優(yōu)位置為Pg=(pg1,pg2,…,pgD),粒子j的速度和位置更新為
其中,ω為慣性權(quán)重,取值范圍為(0,1);c1,c2為學(xué)習(xí)因子,通常c1=c2=2;r1,r2為隨機數(shù),取值范圍為(0,1);粒子位置更新上限為x∈[xmin,xmax],速度更新上限為v∈[vmin,vmax]。
假設(shè)在一個D 維搜索空間中,有Z×n個粒子組成一個多粒子群落,牙齒Ti的粒子群中當前粒子ij的位置為Xij=(xij1,xij2,…,xijD),當前速度為Vij=(vij1,vij2,…,vijD),粒子ij經(jīng)歷的最優(yōu)位置為Pij=(pij1,pij2,…,pijD),粒子群迄今為止經(jīng)歷的最優(yōu)位置為Pig=(pig1,pig2,…,pigD)。
2.2.1 慣性參數(shù)改進
口腔正畸學(xué)中根據(jù)牙齒的功能和特征將牙齒分為:切牙、尖牙及磨牙。從圖4 中可以看出,以牙中線為基準切牙包括中切牙、側(cè)切牙;尖牙包括尖牙、第一雙尖牙和第二雙尖牙,磨牙分為第一磨牙和第二磨牙,根據(jù)牙齒的特性,不同牙齒移動的難易程度不同。臨床表明:在牙齒矯治過程,由于磨牙的生理組織重建困難,其偏移量一般較小,而切牙及尖牙在正畸過程中的偏移量相對會大。
圖4 牙齒分類示意圖 Fig.4 Schematic diagram of tooth classification
在改進多粒子群算法時一個粒子群代表單顆牙齒,各粒子的目標趨于一致,因此可以將ω相應(yīng)增大達到快速收斂、局部最優(yōu)的效果。針對不同牙齒的移動空間和生理結(jié)構(gòu),切牙移動速度最快,其次為尖牙,最后為磨牙。設(shè)置不同的ω值來控制相應(yīng)牙齒的移動速度,使整口牙逐步向理想位置移動,減少了個別牙齒因兩側(cè)牙齒已達到理想位置產(chǎn)生碰撞導(dǎo)致路徑過長的情況,從而間接達到全局優(yōu)化的效果。本文以牙弓深度作為改進的依據(jù),對不同牙齒在粒子群規(guī)劃過程中的慣性權(quán)值ω做個性化設(shè)置。在確定理想位姿時,采用Beta 曲線為標準的牙弓線,牙齒Ti在y軸的值為hi,如圖5 所示。
圖5 牙弓深度示意圖 Fig.5 Schematic diagram of dental arch depth
慣性參數(shù)計算步驟:
步驟1.讀取牙齒STL 數(shù)據(jù),并創(chuàng)建OBB 及整口牙齒中心坐標系;
步驟2.坐標轉(zhuǎn)換,將原始坐標轉(zhuǎn)換到以牙齒中線為y軸,垂直牙平面為z軸,沿后磨牙中心連線為x軸的新坐標系中;
步驟3.根據(jù)牙齒特征點擬合Bate 曲線,擬合理想牙弓線;
步驟4.根據(jù)牙弓線計算牙齒理想位置,根據(jù)牙齒中心位置計算牙齒hi;
步驟5.計算牙齒Ti的慣性參數(shù)ωi用于Mutil-PSO 算法中的速度更新公式。
對不同牙齒設(shè)置不同的慣性權(quán)值,改進后的ωi為牙齒Ti的慣性參數(shù),即
其中,h為牙弓線的最大深度,即Beta 曲線y軸的最大值;hi為第i顆牙齒中心點y軸的值,因此改進后牙齒Ti粒子群的粒子速度和位置更新式為
2.2.2 變步長改進
牙齒在不同的正畸階段,其移動距離在0.1~0.5 mm 之間,付敬鼎[16]在RRT 算法及其他牙齒正畸路徑規(guī)劃中,采用最大值0.5 mm 作為每個階段的移動路徑的搜索范圍,有利于算法快速收斂,但搜索范圍過大會造成頻繁碰撞和路徑曲折過長。為了減少碰撞提高算法的效率,本文通過變步長的方法改進正畸路徑規(guī)劃過程,對于位置更新的初始值設(shè)置為最大值0.5 mm,在其滿足終止條件時則修改步長:
設(shè)置初始步長Sim=0.5 mm,當F1≤n×0.3 mm時,調(diào)整約束Sim=0.3 mm,即設(shè)置位置更新的上限為0.3 mm;當F1≤n×0.1 mm 時,調(diào)整約束Sim=0.1 mm,即上限為0.1 mm;當F1≤n×0.1 mm×0.2 時,終止規(guī)劃。
2.2.3 適應(yīng)度函數(shù)構(gòu)造
在牙齒正畸路徑規(guī)劃中,根據(jù)優(yōu)化目標為最短路徑和最小旋轉(zhuǎn)量,約束條件為在規(guī)劃時不產(chǎn)生碰撞,可構(gòu)造出相應(yīng)的適應(yīng)度函數(shù)
其中,第1 項為碰撞的罰函數(shù);第2 項為牙齒位移量;第3 項為牙齒旋轉(zhuǎn)量;λ1,λ2,λ3為相應(yīng)的權(quán)重,本文分別取100,6,4。
2.2.4 改進算法流程
患者牙齒個數(shù)為n,矯治階段為m,每顆牙齒以一個單獨的粒子群算法進行規(guī)劃,因此每顆牙齒的粒子的編碼為xiyiziαiβiγi,改進多粒子群算偽代碼見算法1。
偽代碼說明:首先對每個粒子群中的粒子進行編碼,初始化種群并設(shè)定相關(guān)參數(shù),其中包括患者牙齒個數(shù)n,需要矯治階段數(shù)m,種群大小Z,每個粒子維數(shù)D,最大迭代次數(shù)M;其次輸入每顆牙齒的初始位置和理想位置,計算各個粒子群的慣性參數(shù),進行路徑規(guī)劃;最后每階段結(jié)束后計算是否達到步長修改條件或終止條件,未達到則繼續(xù)規(guī)劃,否則修改相應(yīng)參數(shù)或終止規(guī)劃。
算法1.多粒子群算法流程。
輸入:牙齒初始初始位姿和牙齒理想位姿。
輸出:牙齒正畸路徑。
該實驗開發(fā)硬件環(huán)境為:Intel i7 1.80 GHz CPU,8 G 內(nèi)存,軟件環(huán)境為:Microsoft Windows 10 操作系統(tǒng)、隱形矯治系統(tǒng)Orthodontic,Matlab開發(fā)環(huán)境和VTK 工具包。
該實驗涉及到擬合Beta 和Spee 曲線獲取牙齒理想位置、對已分割的牙齒建立數(shù)學(xué)模型、牙齒OBB 建立及分離軸碰撞檢測、牙齒運動路徑規(guī)劃方法的實現(xiàn)以及優(yōu)化等過程。本文共對10 口牙齒數(shù)據(jù)進行實驗,選取其中3 組下頜牙齒正畸前和正畸后的位姿數(shù)據(jù),在VTK 中對牙齒位置進行模擬展示如圖6 所示,左側(cè)和右側(cè)分別為正畸前、后的效果圖。經(jīng)對比可以看出,正畸后牙齒排列有明顯的變化,且牙列整齊、美觀,符合正畸的要求。
圖6 牙齒正畸前后對比 Fig.6 Comparison of teeth before and after orthodontic treatment
本文以一組下頜牙齒正畸過程中的6 個階段為例,展示如何獲取牙齒正畸過程中的關(guān)鍵中間位置。圖7 為改進方法的矯治路徑三維可視化效果,其中紅色為牙齒理想位置的OBB,黑色為當前階段牙齒OBB 位置。通過實驗,可以看出正畸的每一階段牙齒都會根據(jù)牙齒運動量約束朝著目標位姿運動,直至達到算法終止條件,該方法能夠為牙齒運動規(guī)劃出無碰撞最優(yōu)路徑,生成符合臨床正畸要求的隱形矯治方案。
圖7 牙齒正畸的6 個階段 Fig.7 Six stages of orthodontic treatment
3.2.1 實驗效率對比
文獻[17]選用多目標粒子群算法,其中一個粒子代表所有牙齒全階段的矯治結(jié)果,即每個粒子由所有牙齒的平移分量和旋轉(zhuǎn)分量組成,因此會導(dǎo)致維度過高;采用粒子運動的位置更新上下限采用生理所能承受最大距離,不利于后期快速朝理想位置靠近。本文采用多粒子群算法進行規(guī)劃,且一個粒子代表一顆牙齒的旋轉(zhuǎn)分量和平移分量,降低了算法的維度,在不同的正畸階段設(shè)置新的位置更新上下限,可以有效縮小搜索范圍,并加快向理想位置靠近。對比實驗中,參數(shù)維度D=6,最大迭代次數(shù)M=50,種群粒子數(shù)Z=50;文獻[17]PSO 算法中ω=0.8,c1=c2=2,r1,r2為隨機數(shù);改進Mutil-PSO算法中,ωi=|hi/h|,c1=c2=2,r1,r2為隨機數(shù),對比實驗中2 算法的終止條件均為F1≤n×0.1 mm×0.2。實驗結(jié)果對10 口牙齒下頜正路徑規(guī)劃的效率求平均值:采用文獻[17]PSO 的粒子群算法,其平均規(guī)劃時間為147.496 403 s,本文算法的平均規(guī)劃時間為91.676 100 s,因此本文方法在獲取牙齒正畸過程中關(guān)鍵中間位置同時也提高了程序運行效率約38%。
3.2.2 正畸效果對比
在牙齒正畸過程中,牙齒從初始位置移動到目標位置移動量和旋轉(zhuǎn)量的大小直接影響患者在正畸過程中所承受的痛苦,因此移動量和旋轉(zhuǎn)量越小越符合正畸的要求。
若牙齒正畸理想位置的旋轉(zhuǎn)中心為Ci(p)=(Cxi(p),Cyi(p),Czi(p)),正畸之后的旋轉(zhuǎn)中心為Ci(m)=(Cxi(m),Cyi(m),Czi(m)),通過比較理想位置和正畸后位置的誤差E來評價正畸效果
由式(5)可知,n顆牙齒m個階段牙齒的位移動量和旋轉(zhuǎn)角度分別表示為F1和F2,同時對比位移誤差E,分析正畸效果和正畸過程中所承受的痛苦,表1 選取6 組數(shù)據(jù)下頜正畸的結(jié)果作為樣本,將本文方法與文獻[17]算法進行對比,其中:位移量總和為F1,正畸誤差為E。
表1 正畸效果對比 Table 1 Comparison of orthodontic effect
分析實驗數(shù)據(jù),對于總旋轉(zhuǎn)角度F2,本文方法與對比方法的差值小于2°;此處主要對比正畸E的值,本文方法正畸后的位置更接近理想位置,正畸效果更好;但在6 組下頜14 顆牙齒的正畸模擬中,位移量稍大,但相應(yīng)的正畸效果明顯好于其他效果,改進前后效果如圖8 所示。
圖8 牙齒正畸效果 Fig.8 Comparison of orthodontic effect
針對虛擬口腔正畸治療系統(tǒng)中牙齒移動路徑規(guī)劃問題,本文基于OBB 提出改進慣性參數(shù)和變步長的多粒子群算法用于牙齒正畸路徑規(guī)劃方法,并在Matlab 中進行仿真試驗,實驗結(jié)果表明改進算法比單粒子群算法效率顯著提升,且更貼近理想位姿的正畸效果。但本文的不足是對牙齒存在缺失等特殊情況未做處理,同時OBB 的碰撞檢測仍然需要耗費大量的時間,因此其矯正方案需要進行更深入的研究,進一步完善和優(yōu)化。