基于改進多粒子群的牙齒正畸路徑規劃

馬 天，楊 秦，李占利

(西安科技大學計算機科學與技術學院，陜西 西安 710054)

近年來，隨著計算機輔助設計仿真、三維重建技術與快速成型技術的快速發展，無托槽隱形矯治引起了廣大正畸醫生以及患者的格外關注，其具有美觀、舒適、安全、便捷等優點[1]。虛擬牙齒正畸系統實現了在三維可視化狀態下牙齒移動的方式和步驟的設計，主要是為制作隱形矯治器提供數據源[2]。該系統首先通過先進的掃描設備獲取患者牙齒的牙頜數字模型[3]，其次在虛擬牙齒正畸系統中進行單顆牙齒分割修補[4]、牙齒移動和矯治方案制定、牙齦組織變形等[5]環節，最終輸出隱形矯治器的母版。其中牙齒移動和矯治方案的確定是非常重要的技術環節，通過計算機算法實現正畸中間路徑的獲取，可以減少對口腔醫師經驗的依賴，提高正畸的精準度，因此牙齒矯治自動化是口腔醫學發展的趨勢[6]。

牙齒正畸是三維環境下有障礙物的多目標路徑規劃問題[7]，傳統的路徑規劃優化方法在應對復雜環境時會存在一定的缺陷[8]。在處理復雜環境信息路徑規劃時，來自于自然界的啟示往往能起到很好的作用，常用的仿生學算法有：蟻群算法、神經網絡算法、粒子群算法和遺傳算法等，其中粒子群采用實數編碼，更方便問題的映射。

國內外有關自動牙齒正畸路徑規劃的研究成果較少，MOTOHASHI 和KURODA[9]的方法首先計算每顆牙齒的特征線段，通過牙弓線牽引的方式實現自動排牙技術，該方法不適宜牙列擁擠碰撞的情況，因此在實際正畸中存在諸多問題；LIN[10]提出半自動排牙技術，仍需提取單顆牙齒特征線段，利用二次函數建立擬合牙弓線，未規避牙齒間可能產生的碰撞和牽引移動過程中牙齒的生理特征，采用該正畸方案可能會帶來負面效果。從1998 年開始，國外開始出現隱形矯治器虛擬牙齒正畸系統，ALIGN 公司[11]自主研發的Invisalign 在正畸臨床上取得了巨大成功；1999 年，GeoDigm 公司[12]開發出一款可以操作牙齒分割、移動及旋轉的虛擬牙齒正畸軟件Emodel；2003 年，美國Cadent 公司[13]實現了牙齒正畸參數分析和牙齒正畸自動排牙算法。國外關于虛擬牙齒正畸系統的文獻大多為公司相關產品的宣傳文案，而有關虛擬牙齒正畸系統技術層面的文獻很少。

近年來，國內在牙齒正畸方面也取得了一些成果，2010 年LI 等[14]采取遺傳算法求解牙齒正畸路徑規劃，認為遺傳算法支持多個物體同時移動符合牙齒正畸的需求，但遺傳算法存在多參數的自適應性問題且效率較低；2016 年山東大學張筱[15]提出單顆牙齒的排牙方法和牙齒正畸過程中的碰撞檢測方法，并設計牙齒正畸路徑規劃的整體思路流程；2018 年付敬鼎[16]提出基于改進快速搜索隨機樹(rapidly-exploring random trees，RRT)算法進行路徑規劃，由于算法本身的隨機性，生成的路徑比較曲折，甚至出現繞遠路。2020 年徐曉強等[17]采用粒子群的算法，以一個粒子代表整口牙齒進行矯治，存在維度過高導致運算效率較低的問題。

根據臨床經驗，牙齒正畸需要考慮其生理組織結構，在以往的研究中，未對不同種類的牙齒進行區分，導致實際和理想正畸過程存在一定的差異。本文改進多粒子群牙齒正畸路徑規劃的方法考慮了多顆牙齒正畸的高維度和碰撞問題，將每顆牙齒采用一個粒子群的方式進行矯治來降低維度以提高算法的效率；并考慮到牙齒的生理組織結構，改進不同粒子群的慣性權值和不同正畸階段位置更新的上下限，也達到了減少碰撞和優化正畸效果的目的。

1 牙齒正畸問題建模

無托槽隱形矯治需要解決牙齒理想位姿確定、牙齒正畸路徑規劃、牙齦組織變形等諸多技術難題，涉及計算機圖形學領域的諸多分支，是一個多學科相結合的復雜系統。本課題主要研究的是在牙齒初始位置和目標位置均已確定的情況下，獲取牙齒正畸中間階段的位置。因此圍繞牙齒矯治過程，將解決實際牙齒三維環境下有障礙的路徑規劃為本文的研究重點。

1.1 牙齒目標位置確定

牙齒正畸路徑規劃是獲取牙齒從初始位置到理想位置的過程，因此確定理想位置的排牙技術是無托槽隱形矯治技術的首要前提。牙齒排列是指在牙齒初始位置已知的情況下，參照排牙原則及理想的牙弓曲線[18]，經過排牙處理得到牙齒的最終位置。

牙齒的理想位置是矯治的重要參考，正確地確定患者牙弓形態極為重要，其是建立穩定、緊密咬合關系的重要基礎。許多專家研究了人類牙弓線數學模型，觀點不一，主要用到的數學模型有橢圓線函數、拋物線函數、多項式函數、Beta方程等[19]。本文選用Beta方程作為擬合函數模型，由上下頜牙齒特征點集在橫斷面投影點擬合曲線，計算各投影點到曲線的最短距離得到平移量Δx，Δy；平移量Δz，則是將牙齒特征點集投影至矢狀面，并通過擬合計算得到，其反映的Spee 曲線[20]的走勢；最終通過計算的平移量來確定理想位置，圖1 為牙齒排列的曲線擬合過程。

圖1 牙齒排列技術((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)

牙齒的旋轉值在口腔醫學中被稱為轉矩角、軸傾角、旋轉角[21]，通過牙齒在空間中任意姿態的變換可由歐拉角構造旋轉變換矩陣求出3 個旋轉量，圖2 為牙齒旋轉示意圖。

圖2 牙齒旋轉分量示意圖((a)歐拉角坐標；(b)牙齒旋轉坐標) Fig.2 Schematic diagram of tooth rotation component ((a) Euler angular coordinates;(b) Spatial coordinates of tooth rotation)

1.2 牙齒碰撞檢測問題分析

牙齒碰撞是指多顆牙齒在正畸路徑上發生交叉時可能出現的情形，牙齒正畸治療的一大難點是如何避免牙齒間的碰撞。由于牙齒本身的不規則性，且其使用的都是大數據量的三角網格，碰撞檢測需要消耗大量時間，因此需要對正畸方案中的碰撞檢測進行調整和簡化。

基于包圍盒的碰撞檢測算法是一類重要的檢測算法[22]，其原理是在三維圖像場景中利用簡單的幾何體即包圍盒包圍碰撞檢測的幾何對象，然后將包圍盒進行碰撞檢測的結果用以替代幾何對象的碰撞檢測結果。包圍盒的形狀越簡單，碰撞檢測算法的效率越高，若包圍盒形狀越接近于研究幾何對象，其越接近真實幾何對象的碰撞檢測結果，所以應盡量選擇形狀簡單且接近幾何對象的包圍盒。常見的包圍盒算法有軸對齊包圍盒(axis-aligned bounding box，AABB)、包圍球、方向包圍盒(oriented bounding box，OBB)以及固定方向凸包(fixed directions hulls，FDH)[23]。牙齒正畸由于空間較為狹窄，需選用盡可能緊密的包圍盒十分必要，且OBB 相交檢測時只需要檢測15 條分離軸，算法復雜度低，因此本文采用OBB 進行碰撞檢測，圖3為采用MATLAB 創建的OBB。

圖3 OBB 創建圖 Fig.3 OBB creation diagram

1.3 牙齒路徑規劃約束條件及目標

對牙齒進行矯治時，需要考慮如生理組織結構重建、牙齒能承受的正畸力范圍、個體因素等，由于這些條件的限制，牙齒在正畸時需要分階段進行，并且每個正畸階段的移動量有限。在理想的正畸階段中，各顆牙齒移動時應該感受不到相鄰牙齒的存在，各自朝著預定的路徑移動，不會發生碰撞。

設牙齒為Ti(i=1,2,…,n)，患者牙齒個數為n，矯治階段為m，牙齒正畸是牙齒從初始姿態Pi0到目標姿態Pim的碰撞路徑，牙齒在第m階段的旋轉中心為Ci(m)=(Cxi(m),Cyi(m),Czi(m))，姿態角為δim(m)=(αi(m),βi(m),γi(m))，牙齒Ti的第m個正畸階段位姿可以表示為Pim=(Ci(m),δi(m))，則牙齒Ti第m個正畸階段的路徑規劃長度為

為了簡化牙齒旋轉的約束條件,則牙齒Ti第m個正畸階段旋轉角度為

根據臨床正畸經驗，牙齒在正畸的一個階段中產生的平移或旋轉的角度不能超過一定范圍，單階段牙齒平移的值Sim不能超過0.5，δim不能超過3°。

牙齒在正畸過程中，若正畸方案不當，可能使牙齒與其相鄰牙齒發生碰撞干涉，影響正畸效果。因此考慮碰撞約束，能夠確保牙齒沿正畸路徑移動時不與相鄰牙齒發生碰撞干涉，設Kim為判斷牙齒Ti在第m階段是否發生碰撞

為了讓相鄰牙齒發生碰撞時，碰撞罰函數Fim能夠有明顯的變化，則給定一個較大常數項L，發生碰撞時罰函數為

在牙齒正畸時，牙齒之間間隙很小，相應的活動空間也很小，因此牙齒正畸路徑規劃問題的優化指標有：所有牙齒移動路徑最短、旋轉角度最小等；且需滿足式(4)的碰撞約束。由式(1)和(2)可得，n顆牙齒m個階段牙齒的最短移動路徑和最小旋轉角度分別為

2 基于改進粒子群的正畸路徑規劃

2.1 粒子群算法

粒子群算法初始化為一群隨機的粒子(隨機解)，根據迭代找到最優解。所有的粒子都有一個由被優化的函數決定的適應值，每個粒子還有一個速度決定其飛出的方向和距離，然后粒子就追隨當前的最優粒子在解空間中搜索。粒子具有速度和位置2 個屬性，速度代表移動的快慢，位置代表移動的方向。每個粒子單獨搜尋的最優解叫做個體極值，粒子群中最優的個體極值作為當前全局最優解，不斷迭代更新速度和位置，最終得到滿足終止條件的最優解。

假設在一個D 維搜索空間中，有Z個粒子組成一個群落，其中第j個粒子表示為一個D 維的向量，則粒子j當前位置為Xj=(xj1,xj2,…,xjD)，當前速度為Vj=(vj1,vj2,…,vjD)，該粒子j經歷的最優位置為Pj=(pj1,pj2,…,pjD)，該粒子群的最優位置為Pg=(pg1,pg2,…,pgD)，粒子j的速度和位置更新為

其中，ω為慣性權重，取值范圍為(0,1)；c1，c2為學習因子，通常c1=c2=2；r1，r2為隨機數，取值范圍為(0,1)；粒子位置更新上限為x∈[xmin,xmax]，速度更新上限為v∈[vmin,vmax]。

2.2 算法改進

假設在一個D 維搜索空間中，有Z×n個粒子組成一個多粒子群落，牙齒Ti的粒子群中當前粒子ij的位置為Xij=(xij1,xij2,…,xijD)，當前速度為Vij=(vij1,vij2,…,vijD)，粒子ij經歷的最優位置為Pij=(pij1,pij2,…,pijD)，粒子群迄今為止經歷的最優位置為Pig=(pig1,pig2,…,pigD)。

2.2.1 慣性參數改進

口腔正畸學中根據牙齒的功能和特征將牙齒分為：切牙、尖牙及磨牙。從圖4 中可以看出，以牙中線為基準切牙包括中切牙、側切牙；尖牙包括尖牙、第一雙尖牙和第二雙尖牙，磨牙分為第一磨牙和第二磨牙，根據牙齒的特性，不同牙齒移動的難易程度不同。臨床表明：在牙齒矯治過程，由于磨牙的生理組織重建困難，其偏移量一般較小，而切牙及尖牙在正畸過程中的偏移量相對會大。

圖4 牙齒分類示意圖 Fig.4 Schematic diagram of tooth classification

在改進多粒子群算法時一個粒子群代表單顆牙齒，各粒子的目標趨于一致，因此可以將ω相應增大達到快速收斂、局部最優的效果。針對不同牙齒的移動空間和生理結構，切牙移動速度最快，其次為尖牙，最后為磨牙。設置不同的ω值來控制相應牙齒的移動速度，使整口牙逐步向理想位置移動，減少了個別牙齒因兩側牙齒已達到理想位置產生碰撞導致路徑過長的情況，從而間接達到全局優化的效果。本文以牙弓深度作為改進的依據，對不同牙齒在粒子群規劃過程中的慣性權值ω做個性化設置。在確定理想位姿時，采用Beta 曲線為標準的牙弓線，牙齒Ti在y軸的值為hi，如圖5 所示。

圖5 牙弓深度示意圖 Fig.5 Schematic diagram of dental arch depth

慣性參數計算步驟：

步驟1.讀取牙齒STL 數據，并創建OBB 及整口牙齒中心坐標系；

步驟2.坐標轉換，將原始坐標轉換到以牙齒中線為y軸，垂直牙平面為z軸，沿后磨牙中心連線為x軸的新坐標系中；

步驟3.根據牙齒特征點擬合Bate 曲線，擬合理想牙弓線；

步驟4.根據牙弓線計算牙齒理想位置，根據牙齒中心位置計算牙齒hi；

步驟5.計算牙齒Ti的慣性參數ωi用于Mutil-PSO 算法中的速度更新公式。

對不同牙齒設置不同的慣性權值，改進后的ωi為牙齒Ti的慣性參數，即

其中，h為牙弓線的最大深度，即Beta 曲線y軸的最大值；hi為第i顆牙齒中心點y軸的值，因此改進后牙齒Ti粒子群的粒子速度和位置更新式為

2.2.2 變步長改進

牙齒在不同的正畸階段，其移動距離在0.1～0.5 mm 之間，付敬鼎[16]在RRT 算法及其他牙齒正畸路徑規劃中，采用最大值0.5 mm 作為每個階段的移動路徑的搜索范圍，有利于算法快速收斂，但搜索范圍過大會造成頻繁碰撞和路徑曲折過長。為了減少碰撞提高算法的效率，本文通過變步長的方法改進正畸路徑規劃過程，對于位置更新的初始值設置為最大值0.5 mm，在其滿足終止條件時則修改步長：

設置初始步長Sim=0.5 mm，當F1≤n×0.3 mm時，調整約束Sim=0.3 mm，即設置位置更新的上限為0.3 mm；當F1≤n×0.1 mm 時，調整約束Sim=0.1 mm，即上限為0.1 mm；當F1≤n×0.1 mm×0.2 時，終止規劃。

2.2.3 適應度函數構造

在牙齒正畸路徑規劃中，根據優化目標為最短路徑和最小旋轉量，約束條件為在規劃時不產生碰撞，可構造出相應的適應度函數

其中，第1 項為碰撞的罰函數；第2 項為牙齒位移量；第3 項為牙齒旋轉量；λ1，λ2，λ3為相應的權重，本文分別取100，6，4。

2.2.4 改進算法流程

患者牙齒個數為n，矯治階段為m，每顆牙齒以一個單獨的粒子群算法進行規劃，因此每顆牙齒的粒子的編碼為xiyiziαiβiγi，改進多粒子群算偽代碼見算法1。

偽代碼說明：首先對每個粒子群中的粒子進行編碼，初始化種群并設定相關參數，其中包括患者牙齒個數n，需要矯治階段數m，種群大小Z，每個粒子維數D，最大迭代次數M；其次輸入每顆牙齒的初始位置和理想位置，計算各個粒子群的慣性參數，進行路徑規劃；最后每階段結束后計算是否達到步長修改條件或終止條件，未達到則繼續規劃，否則修改相應參數或終止規劃。

算法1.多粒子群算法流程。

輸入：牙齒初始初始位姿和牙齒理想位姿。

輸出：牙齒正畸路徑。

3 實驗結果及分析

該實驗開發硬件環境為：Intel i7 1.80 GHz CPU，8 G 內存，軟件環境為：Microsoft Windows 10 操作系統、隱形矯治系統Orthodontic，Matlab開發環境和VTK 工具包。

3.1 排牙效果分析

該實驗涉及到擬合Beta 和Spee 曲線獲取牙齒理想位置、對已分割的牙齒建立數學模型、牙齒OBB 建立及分離軸碰撞檢測、牙齒運動路徑規劃方法的實現以及優化等過程。本文共對10 口牙齒數據進行實驗，選取其中3 組下頜牙齒正畸前和正畸后的位姿數據，在VTK 中對牙齒位置進行模擬展示如圖6 所示，左側和右側分別為正畸前、后的效果圖。經對比可以看出，正畸后牙齒排列有明顯的變化，且牙列整齊、美觀，符合正畸的要求。

圖6 牙齒正畸前后對比 Fig.6 Comparison of teeth before and after orthodontic treatment

本文以一組下頜牙齒正畸過程中的6 個階段為例，展示如何獲取牙齒正畸過程中的關鍵中間位置。圖7 為改進方法的矯治路徑三維可視化效果，其中紅色為牙齒理想位置的OBB，黑色為當前階段牙齒OBB 位置。通過實驗，可以看出正畸的每一階段牙齒都會根據牙齒運動量約束朝著目標位姿運動，直至達到算法終止條件，該方法能夠為牙齒運動規劃出無碰撞最優路徑，生成符合臨床正畸要求的隱形矯治方案。

圖7 牙齒正畸的6 個階段 Fig.7 Six stages of orthodontic treatment

3.2 實驗效果對比

3.2.1 實驗效率對比

文獻[17]選用多目標粒子群算法，其中一個粒子代表所有牙齒全階段的矯治結果，即每個粒子由所有牙齒的平移分量和旋轉分量組成，因此會導致維度過高；采用粒子運動的位置更新上下限采用生理所能承受最大距離，不利于后期快速朝理想位置靠近。本文采用多粒子群算法進行規劃，且一個粒子代表一顆牙齒的旋轉分量和平移分量，降低了算法的維度，在不同的正畸階段設置新的位置更新上下限，可以有效縮小搜索范圍，并加快向理想位置靠近。對比實驗中，參數維度D=6，最大迭代次數M=50，種群粒子數Z=50；文獻[17]PSO 算法中ω=0.8，c1=c2=2，r1，r2為隨機數；改進Mutil-PSO算法中，ωi=|hi/h|，c1=c2=2，r1，r2為隨機數，對比實驗中2 算法的終止條件均為F1≤n×0.1 mm×0.2。實驗結果對10 口牙齒下頜正路徑規劃的效率求平均值：采用文獻[17]PSO 的粒子群算法，其平均規劃時間為147.496 403 s，本文算法的平均規劃時間為91.676 100 s，因此本文方法在獲取牙齒正畸過程中關鍵中間位置同時也提高了程序運行效率約38%。

3.2.2 正畸效果對比

在牙齒正畸過程中，牙齒從初始位置移動到目標位置移動量和旋轉量的大小直接影響患者在正畸過程中所承受的痛苦，因此移動量和旋轉量越小越符合正畸的要求。

若牙齒正畸理想位置的旋轉中心為Ci(p)=(Cxi(p),Cyi(p),Czi(p))，正畸之后的旋轉中心為Ci(m)=(Cxi(m),Cyi(m),Czi(m))，通過比較理想位置和正畸后位置的誤差E來評價正畸效果

由式(5)可知，n顆牙齒m個階段牙齒的位移動量和旋轉角度分別表示為F1和F2，同時對比位移誤差E，分析正畸效果和正畸過程中所承受的痛苦，表1 選取6 組數據下頜正畸的結果作為樣本，將本文方法與文獻[17]算法進行對比，其中：位移量總和為F1，正畸誤差為E。

表1 正畸效果對比 Table 1 Comparison of orthodontic effect

分析實驗數據，對于總旋轉角度F2，本文方法與對比方法的差值小于2°；此處主要對比正畸E的值，本文方法正畸后的位置更接近理想位置，正畸效果更好；但在6 組下頜14 顆牙齒的正畸模擬中，位移量稍大，但相應的正畸效果明顯好于其他效果，改進前后效果如圖8 所示。

圖8 牙齒正畸效果 Fig.8 Comparison of orthodontic effect

4 結束語

針對虛擬口腔正畸治療系統中牙齒移動路徑規劃問題，本文基于OBB 提出改進慣性參數和變步長的多粒子群算法用于牙齒正畸路徑規劃方法，并在Matlab 中進行仿真試驗，實驗結果表明改進算法比單粒子群算法效率顯著提升，且更貼近理想位姿的正畸效果。但本文的不足是對牙齒存在缺失等特殊情況未做處理，同時OBB 的碰撞檢測仍然需要耗費大量的時間，因此其矯正方案需要進行更深入的研究，進一步完善和優化。