基于改进多粒子群的牙齿正畸路径规划

马 天，杨 秦，李占利

(西安科技大学计算机科学与技术学院，陕西 西安 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 的碰撞检测仍然需要耗费大量的时间，因此其矫正方案需要进行更深入的研究，进一步完善和优化。