An orthodontic path planning method based on improved gray wolf optimization algorithm

Automatic tooth arrangement and path planning play an essential role in computer-aided orthodontic treatment. However, state-of-the-art methods have some shortcomings: low efficiency, excessive cost of displacements or collisions and insufficient accuracy. To address these issues, this paper proposes an innovative orthodontic path planning method based on the improved gray wolf optimization algorithm, which is called OPP-IGWO. First, the tooth model is preprocessed to obtain initial pose in which each tooth is segmented and built up with an oriented bounding box. Next, the target pose of each tooth is determined through the ideal dental arch curve and the optimal jaw principle. Finally, the path from the initial pose to the target pose of each tooth is planned based on IGWO, which is improved mainly from three aspects: (1) The greedy idea is adopted to initialize the gray wolf population based on dental interpolation. (2) The linear convergence factor in the traditional gray wolf optimization algorithm (GWO) is replaced with a nonlinear convergence factor. (3) We propose a position update strategy based on a dynamic weighting approach, which introduces a learning rate for each gray wolf. The experimental results show that our OPP-IGWO method outperforms the state-of-the-art methods. Compared with improved genetic algorithm, multiparticle improved swarm optimization, normal simplified mean particle swarm optimization algorithm, improved artificial bee colony algorithm, chaotic gray wolf optimization algorithm, hybrid gray wolf optimization algorithm and differential evolution algorithm, random walk gray wolf Optimization algorithm and reinforcement learning gray wolf optimization algorithm, the OPP-IGWO has an improvement on performance by 15.46%, 2.24%, 7.18%, 15.99%, 7.67%, 1.01%, 1.42%, 10.23%, respectively.


Introduction
Nowadays, with the rapid development of stomatology in the world, orthodontic automation is becoming the trend of this industry, such as automatic 3D dental scan mapping (Woodsend et al. 2021), automatic generation of the spinal canal contour of the tooth (Huo et al. 2020), automatic prediction of the visual effect of orthodontic treatment in portrait images (Chen et al. 2022) and so on.Under this trend, the automation of orthodontic path planning realized by computer technology can not only quickly help orthodontists predict the final orthodontic effect and assist in the design of orthodontic schemes, but make the patients know the whole process of orthodontic treatment at a glance through the generated 3D animation scheme.Therefore, orthodontic path planning is the core technology of invisible orthodontic scheme generation.The research on the automation of orthodontic path planning has great clinical value.
At present, most dental technicians still use the more complex human-computer interactive planning method (Motohashi and Kuroda 1999;Nota et al. 2022), which usually takes an average of 30 min to complete a patient's orthodontic plan design.Over the recent years, some studies on orthodontic path planning have been taken, and these methods have addressed some issues.However, they still have some shortcomings to be improved.The researches of orthodontic path planning are mainly divided into two directions.One direction is to learn from the path planning of robot movement.These methods use more traditional heuristic search algorithms, such as A*, Dijkstra, and RRT (Li et al. 2019) algorithms to search for tooth paths.However, due to the strong randomness of the algorithms, this kind of methods may make the generated tooth path too tortuous, or even take a long detour.Another direction is to improve the classical bionic intelligence algorithms and apply them to the field of tooth path planning, such as GA (Li et al. 2009), PSO (Xu et al. 2020), ABC (Li et al. 2020) and so on.However, some of these methods tend to fall into local optimality when faced with high-dimensional problems such as orthodontic path planning, while others do not adequately consider the physiological conditions in the mouth and generate unrealistic orthodontic paths.So, the effect of generating tooth paths is not ideal.
To address these issues, this paper proposes a new orthodontic path planning method based on the improved gray wolf optimization algorithm, called OPP-IGWO.First, the tooth model is preprocessed to obtain initial pose in which each tooth is segmented and built up with an oriented bounding box (OBB).Then, through the ideal dental arch curve and the optimal jaw principle, the target pose for orthodontic treatment is determined.Finally, our IGWO method is used to plan the orthodontic path of the teeth and generate an invisible orthodontic plan that meets the clinical requirements.
The main contributions of this paper are as follows: 1.An improved grey wolf optimization algorithm (IGWO) is proposed and applied to the field of orthodontic path planning for the first time, which is more effectively to obtain the optimal orthodontic path.The traditional GWO algorithm is improved mainly from three aspects: (1) In order to avoid each gray wolf searching inefficiently in the initial stage of the search space, the greedy idea is adopted to initialize the gray wolf population based on dental interpolation.And the global optimal solution is split into m local optimal solutions, which improves the search ability of the algorithm in the initial phase.(2) The linear convergence factor in the traditional GWO is replaced with a nonlinear convergence factor, which improves the global search and local mining ability.(3) We propose a position update strategy based on a dynamic weighting approach, which introduces a learning rate for each gray wolf to improve the global search ability.
2. Beta curve and Spee curve are introduced to more accurately fit the ideal dental arch curve.At the same time, the offset of a single tooth in the 6D degree of freedom is determined by referring to the principle of optimal jaw.Then, the target pose of each tooth is determined through the ideal dental arch curve and the optimal jaw principle, which reduces the time consumption.
The rest of this paper is organized as follows.Section 2 reviews related work and discusses the unresolved issues.In Sect.3, the proposed OPP-IGWO method is described in detail.In Sect.4, some experiments are conducted to evaluate the feasibility and effectiveness of our method.Section 5 concludes this paper and discusses some future work.

Related work
In recent years, scholars have carried out researches on the automation of generating orthodontic plans in invisible orthodontic systems and proposed some feasible methods.For the automation of orthodontic path planning, the target pose of each tooth is generated first, and then, the path from the initial pose to the target pose of each tooth is simulated by computer.Therefore, the target pose generation methods are first discussed, and then, the orthodontic path planning methods are described.

Target pose generation methods
According to the original model of the patient's teeth and the principle of orthodontic arrangement, the final state of the dental model after orthodontic treatment is preferentially obtained through the operation of translation or rotation.We call this process the determination of the tooth target pose.The commonly used method is manual interactive tooth arrangement, which was first proposed by Motohashi and Kuroda (1999).According to the actual situation of the patient's dental arch and malocclusion, the orthodontist manually arranges teeth by human-computer interaction in the invisible orthodontic software and adjusts the pose of the patient's teeth to the ideal state one by one.Alan (2005) proposed a semi-automatic tooth arrangement method, which defined a characteristic line segment in the direction of the tooth crown for each tooth by fitting the dental arch line with a quadratic function, and guided the tooth to approach the target line segment of the dental arch by using the position of the line segment on the dental arch line.However, this method does not consider the interdental collision, which may cause the patient's condition to deteriorate.Kumar et al. (2013) proposed a spring mass model based on tooth surface features.By defining tooth surface feature constraints such as cusps and edge ridges, the ideal occlusal relationship of the upper and lower jaws was established.However, this method has too many iterations and too much running time.Dai et al. (2018) proposed a reconfigurable rule-driven automatic tooth arrangement method for complete denture.Four typical operators are used to establish the constraint mapping association between teeth and patient constraint sets.Different clinical tooth arrangement rules can be flexibly realized by using the process reorganization of different constraint operators.However, this method needs to adjust the proportion parameters of artificial teeth to get a better shape.Wei et al. (2020) developed the automatic tooth arrangement framework TANet based on deep learning, which effectively solves the problem of structural 6D pose prediction generated in the tooth arrangement target.However, the input tooth model relies too much on the training data.At the same time, the tooth target pose predicted by the neural grid does not take into account the constraints of oral physiology, so it needs further processing.
The method used in this paper to determine the target pose is different from the aforementioned methods.First, the feature points of each tooth are extracted.Then, the Beta curve and the Spee curve are introduced to fit the ideal dental arch curve, and referring to the principle of optimal jaw, the offset of a single tooth in 6D degrees of freedom is calculated to generate the target pose of the tooth model.

Orthodontic path planning methods
After obtaining the target pose for orthodontic treatment, the orthodontist can plan the path from the initial pose to the target pose.Due to the limitations of various factors such as the ability of the physiological tissues of the teeth to bear the orthodontic force, the orthodontist needs to design the orthodontic plan in stages and then uses 3D printing technology to make the tooth models in each stage into a series of invisible appliances (Yu et al. 2022).Li et al. (2009) used the genetic algorithm (GA) to solve the tooth path problem.The method is feasibility, but the genetic algorithm is inefficient and prone to premature maturity.Li et al. (2009) proposed an improved RRT algorithm for tooth path search, which realized the collision-free path search.But the randomness of the algorithm made the generated paths tortuous.Liu et al. (2020) improved the artificial bee colony algorithm (ABC).The method uses the evaluation indicators of the teeth as the objective function and obtains the optimal path of the teeth by using external storage to update the Pareto solution.Ma et al. (2021) considered that the previous tooth path planning did not take into account that different types of teeth have different degrees of correction.So, they improved the particle swarm algorithm by adding the depth ratio of the dental arch to the inertia parameter w to control the search speed.And the method introduced the idea of multi-particle swarm (Mutil-IPSO) and changed the way of updating the upper and lower limits of the position to realize the refined orthodontic treatment of the teeth and the automatic generation of the orthodontic plan.Although the parameters of the particle swarm algorithm are simple and easy to implement, its performance is not good when dealing with high-dimensional problems such as the orthodontic path planning.Moreover, it is easy to fall into the local optimum.
The proposed OPP-IGWO approach in this paper is different from the aforementioned methods.First, the greedy idea is used to initialize the population.The initial position of each gray wolf is generated by tooth interpolation.Next, a nonlinear convergence factor a is introduced as an adjustment strategy to better coordinate the global search and local mining capabilities in the GWO algorithm.Finally, we modify the position update formula of the population.The search step size and direction of gray wolves are adjusted by the learning rate w.

The proposed OPP-IGWO method
The proposed OPP-IGWO method aims to automatically generate an invisible orthodontic plan based on the patient's tooth model, and its procedure is shown in Fig. 1.First, the tooth model needs to be preprocessed to segment individual teeth and establish their respective OBB (You et al. 2022).Then, the ideal dental arch curve is fitted by the Beta curve and the Spee curve to determine the target pose for orthodontic treatment.Finally, the orthodontic path is planned by our IGWO method.

Tooth model preprocessing
In this step, first, we use seed point diffusion, watershed algorithm, etc. techniques to separate the tooth and gingival models in the dental-jaw model.The tooth space information is retained during segmentation, and the side of each tooth is repaired after the segmentation.Then, a single tooth model with a natural shape can be obtained.The tooth models in this paper are numbered using the Fe ´de ´ration Dentaire Internationale (FDI) tooth position marking method (ALShami et al. 2019), which is commonly used worldwide.Each tooth is numbered using 2 Arabic numerals with the tens representing the quadrant of the tooth and the digits representing the position of the tooth.Figure 2 shows a diagram of the FDI tooth position marking method and shows the numbers of the teeth in the four regions of the patient's mouth.
Second, in order to avoid collisions during orthodontic treatment, it is necessary to introduce collision detection between teeth.The tooth models studied in this paper are stored in STL format (Shah et al. 2018) and composed of a large number of triangular meshes.It will take a lot of time to perform collision detection directly, so it is also necessary to simplify the tooth model to improve the detection efficiency.
Due to the rigid body properties of teeth, two teeth cannot interfere with each other or even penetrate each other in the same space and time.Therefore, we establish an oriented bounding box (OBB) (You et al. 2022) for the segmented teeth, respectively, and use the separation axis theorem (Liang et al. 2015) to quickly judge the collision between teeth.For any group of OBB in three-dimensional space, there are 15 potential separation axes, of which 6 are the three-axis directions of bounding boxes A and B, respectively, and the remaining 9 are the cross-product directions of any two coordinate axis vectors in bounding boxes A and B. Therefore, for a set of OBB, only 15 tests are required at most.If there is a separation axis that meets the condition, it can be judged that the two do not collide.Taking a young female patient P as an example, Fig. 3 shows a schematic diagram of the establishment of OBB of the 14 teeth of patient P's lower jaw.

Determination of the target pose of each tooth
Orthodontic path planning is the process of obtaining the path from the initial pose to the target pose of teeth so determination of target pose of each tooth is the primary prerequisite for orthodontic path planning.The current method of automatic tooth arrangement mainly refers to the ideal dental arch curve in the patient's mouth to guide the arrangement of teeth (Zou et al. 2012).In this paper, first, we select the feature points of each tooth.Next, the Beta curve and the Spee curve are used to fit the ideal dental arch curve (Stanley et al. 1998).Finally, the displacement and rotation of each tooth are calculated by the marked feature points to determine the target pose.

Select the feature points of the teeth
In this step, we mark each tooth to record its pose state during movement.Tooth feature points are usually distributed on the surface of the crown, such as the wide points on both sides of the incisor, the apex of the canine, the four prominences of the molar crown and the lowest point of the gum line.We select the feature points of different types of teeth in turn.Figure 4 shows the feature points of four different types of teeth.

Calculate the displacement of the tooth
As shown in Fig. 5, we can split the displacement of dental treatment into horizontal displacement and vertical displacement.
For solving the horizontal displacements Dx and Dy, all the feature points of the upper and lower teeth are projected on the jaw plane (XOY).According to the fitting method of the Beta curve (Stanley et al. 1998;Levent et al. 2019), we select the feature points of the incisal edge of the bilateral central incisors, the cusps of the canines and the distal buccal cusps of the second molars to complete the fitting.The mathematical expression of the Beta curve is shown in Eq. 1: where D is the depth of the dental arch, W is the width of the dental arch and e is the curvature parameter.
Next, we need to calculate the position of the shortest distance of each projection point to the ideal dental arch curve.The distance is between the projection of the feature point and the tangent of the curve.The calculation formula is shown in Eq. 2: where x 0 is the x-coordinate of the above-selected initial feature point projected on the XOY plane.The  For solving the vertical displacement Dz, all the feature points of the upper and lower teeth are projected to the sagittal plane of the oral cavity (YOZ).At this time, the projection trend of the feature points is close to the Spee curve (Vieira et al. 2021;Paes-Souza et al. 2022) of the teeth, so the Spee curve can be fitted by a quartic polynomial, and the curve equation is shown in Eq. 3: Since the Spee curve is located on the sagittal plane of the lateral surface of the oral cavity, Dz can be obtained by substituting the Dy obtained above into Eq.3. In this way, the total displacement of a tooth before and after orthodontic treatment in the tooth row coordinate system is determined.

Calculate the rotation of the tooth
In orthodontics, the rotation of teeth includes torque angle a, axial inclination angle b and torsion angle c, which correspond to the angles generated by the rotation of the teeth around the three axes of X, Y and Z in the local coordinate system, respectively, as shown in Fig. 6.
For the torque angle a and the axial inclination angle b, according to the description in the six criteria of Andrews normal jaw and the optimal dentition criteria proposed by some researchers (Chen et al. 2021), this paper directly selects the optimal inclination angle of the teeth in the ideal dentition, as shown in Table 1.
From Table 1, the torque amount Da and the axial inclination amount Db of a tooth before and after orthodontic treatment can be obtained.For the twist angle c, the Y axis of the local coordinate system of a single tooth in the ideal dentition is perpendicular to the tangent direction of the ideal dental arch curve, and the twist amount Dc can be obtained.In this way, the total rotation amount of a tooth before and after orthodontic treatment in the tooth row coordinate system is determined.
The result of the target pose obtained by performing the above operation on the tooth model of patient P is shown in Fig. 7. From before and after the teeth arrangement of the upper and lower jaw, after algorithm iterative adjustment, all the teeth of patient P are adjusted to the target pose through local displacement and rotation and are neatly and closely arranged according to the ideal dental arch curve.From the perspective of before and after the teeth arrangement of whole jaw, the malocclusion in the anterior teeth area before the teeth arrangement is more serious, but it is greatly improved after the teeth arrangement.

Improved gray wolf optimization algorithm
After the initial pose and target pose of the teeth are obtained, the orthodontist can use the virtual dental orthodontic system to further design an appropriate invisible orthodontic solution for the tooth movement.The core of designing a patient's orthodontic plan is to plan the intermediate stage of the tooth movement path.In order to improve the accuracy of the path plan and the effect of the generated path, this paper adopts the gray wolf optimization algorithm (GWO) (Mirjalili et al. 2014).The basic idea is to use the fitness function to filter out the three head wolves a,b, d, who lead the remaining wolves x to move toward them.The gray wolf population search for the global optimal solution in the continuous population iteration.However, the traditional gray wolf optimization algorithm is easy to fall into local optimum when solving complex optimization problems, which results in low accuracy of results 0 et al. 2021).Therefore, we propose an improved gray wolf algorithm (IGWO), which improves the traditional gray wolf optimization algorithm from three aspects: population interpolation initialization, nonlinear convergence factor and weighted position update.

Limitations analysis
Although the gray wolf optimization algorithm has the advantages of simple formulas, fewer parameters, higher robustness and better ability to jump out of the local optimum, there are still several drawbacks when the GWO is used to solve the orthodontic path planning problem (Liu et al. 2021).

Population random initialization
When the GWO initializes the gray wolves, if there is no relevant constraint, the GWO will lead to randomly generated dental path points in a messy and disorderly manner, at which time the fitness values of the initial gray wolves are too poor to converge easily, which greatly affects the search efficiency.An orthodontic path planning method based on improved gray wolf optimization algorithm 16595 123

Linear convergence factor
The balance between the global and local search capabilities of GWO is dynamically controlled by the stochastic parameters.The convergence factor in the stochastic parameters is linearly decreasing, but orthodontic path optimization is a nonlinear process.The linear convergence factor cannot accurately model the nonlinear optimization process, which results in poor accuracy in the face of complex optimization problems.

Balanced position update
The gray wolf position update is influenced by the average of the three individuals a; b and d with the highest fitness values, which does not reflect the high and low status of the three individuals.And if the initial position of the three individuals is poor, the optimization search process can easily fall into the local optimum and stagnate, which leads to premature convergence.For complex optimization problems, it is generally required that the optimization algorithms can traverse the range of possible global optimal solutions as much as possible, while the position update strategy of the GWO is not flexible, which makes it difficult to search for the global optimal solution when solving some complex problems.

Improvement strategies
We improved the GWO in the following three aspects to make it applicable to orthodontic path planning.

Population interpolation initialization
In order to avoid each gray wolf searching inefficiently in the initial stage of the search space, we introduce the greedy idea to initialize the gray wolf population.Orthogonal path planning is a high-dimensional problem, so this paper uses interpolation methods to reduce the computational effort (Aydınlılar et al. 2021;Vyas et al. 2021).In the orthodontic path planning process, if the constraints such as collision are not considered, the ideal intermediate stage path found in this paper is actually the movement of the teeth directly from the initial pose to the  The interpolation method for initializing gray wolf population based on greedy idea is shown in Eq. 4: where C si and C ei represent the initial pose and target pose of the teeth, respectively.U is the uniform distribution function of the 6-dimensional coordinates containing the displacement and rotation information.The initial solution of X i obeys a uniform distribution in the interval, which improves the initial fitness average of the gray wolf population.

Nonlinear convergence factor
The population intelligence optimization algorithm needs to coordinate the global search ability and the local exploitation ability in the search process.The global search seeks the breadth of search and explores the possible locations of optimal solutions in the solution space as much as possible, which can effectively enhance the diversity of the population and seek the global optimal solution in the early stage of search, while avoiding the premature convergence of the algorithm.The local exploitation seeks the depth of the search and explores the local optimal solution locations as much as possible, which can effectively search for the optimal solution in depth in the later stage of the search.
In the GWO, the parameter a controls the global search and local exploitation ability as shown in Eq. 5, which is not suited to optimize the orthodontic path.Therefore, we improve it by introducing an adjustment strategy with a nonlinear convergence factor a as shown in Eq. 6: where e is the base of natural logarithm, t is the current iteration, T max is the maximum number.The comparison of linear convergence factor and nonlinear convergence factor is shown in Fig. 9.It can be seen from Fig. 9 that when the maximum number of iterations is set to 500, the nonlinear convergence factor a shows a decreasing pattern of fast first and then slow with the increase in the number of iterations.The convergence rate is slower in the early stage of iteration, which is conducive to global search and exploration of global optimal solution; at the end of the iteration, the convergence rate is accelerated, which is conducive to local search, enhances local exploitation ability, and seeks local optimal solutions faster and more accurately.

Weighted position update
In the traditional gray wolf algorithm, the search step size and direction of gray wolves x are determined by the three head wolves a; bandd which are screened out by the fitness  function.But the gray wolves x are guided by the three head wolves a; bandd in a static average, which does not reflect the leadership of the three.Moreover, the position of a is not necessarily the global optimal position.If the three head gray wolves force the gray wolves to keep approaching to them and completely ignore the information exchange between gray wolves x, it is easy to fall into the local optimum and cause the search to stagnate.Therefore, this paper proposes a position update strategy based on a dynamic weighting approach (Song et al. 2022), which introduces learning rate w for gray wolves x as shown in Eq. 7: where w 1 ; w 2 ; w 3 correspond to the learning rate of gray wolves x for heard wolves a, b, d, respectively.X 1 ; X 2 ; X 3 correspond to the current position of heard wolves a, b, d, respectively.It can be seen that the farther the head wolf is from the gray wolves x, the higher the learning rate is, which results in larger search step size of approaching this head wolf, so as to expand the global search ability of the algorithm.
The improved final position update formula of gray wolves x is shown in Eq. 8: where the search step size and direction of the gray wolves are adjusted by the learning rate w i , which improves the ability of gray wolf population to explore the prey globally and attack the prey locally to avoid fall into local optimum due to the extreme circumstances of the head wolves.

Orthodontic path planning
The process of orthodontic path planning based on the IGWO algorithm is shown in Fig. 10, which mainly includes 5 sub-steps.First, the encoding of the position of each gray wolf is determined, which corresponds to the solution of the orthodontic path planning.Then, the initial positions of gray wolf population are generated based on the interpolation initialization strategy of greedy thinking.
Next, the relevant parameters of the IGWO algorithm are initialized according to the improved mathematical formula.Then, we construct a fitness function by objective function and constraints of the tooth path planning problem to calculate the fitness value of the gray wolf population to screen out the optimal three head wolves a, b, d.Then, a position update strategy based on dynamic weighting is used to guide the wolves x to better complete the search.Finally, the optimal orthodontic path is obtained by repeating the above sub-steps.
The details of each sub-step will be described in one by one as follows.

Encode the position
In this step, this paper introduces the coding method of genes in the genetic algorithm to encode the position of each gray wolf.In a gray wolf population, each gray wolf represents a tooth movement arrangement strategy which includes the displacement information ðx ij ; y ij ; z ij Þ and rotation information ða ij ; b ij ; c ij Þ of each tooth at each step.And the search dimension D of the gray wolf population is determined according to the number of variables in the actual problem.If the number of teeth of the patient is n, the number of orthodontic steps is m and each tooth contains six degrees of freedom generated by displacement and rotation at the same time, the search dimension D at this time has a total of 6 9 n 9 m dimensions.Therefore, the gray wolf can be encoded according to the following structure:

Initialize the population
This sub-step is to initialize the gray wolf population.We introduce the greedy idea to split the global optimal solution into m local optimal solutions according to the path points and take the initial solution generated based on dental interpolation as the initial position of the population, which can improve the mean value of the initial fitness of the gray wolf population and avoid the gray wolf blind search to improve the optimization efficiency in the early stage of the search.

Initialize and update related parameters
The third sub-step is to initialize the relevant parameters a, A, and C in the gray wolf optimization algorithm.
The hunting behavior of gray wolf population is defined in Eqs. 9 and 10: where t is the current iteration, D is the distance between the gray wolf and the prey, X p and X are the position vectors of the prey and the gray wolf, respectively, and the expressions of the coefficient vectors A and C are shown in Eqs.11 and 12, respectively: where r 1 ; r 2 2 ½0; 1 are random and a 2 ½0; 2 is the convergence factor.To better coordinate the global search and local exploitation ability in the orthodontic path optimization process, we improve the convergence factor a with nonlinear adjustment as shown in Eq. 5.

Calculate fitness value
The fourth sub-step is to calculate the fitness value of each gray wolf.The construction of fitness function is an important part of the optimization process.On the one hand, the fitness function needs to be closely combined with the actual optimization problem, and the objective function and the constraints in the problem must be mapped in a certain way.On the other hand, the fitness function is also the basis for evaluating the positions of gray wolves in each iteration in the algorithm and used to screen the optimal solution.Orthodontic path planning is to find an optimal path under the condition of orthodontics.According to the disassembly of the motion state, there are two main goals: one is the minimum sum of displacements and another is the minimum sum of rotations.In addition, two main constraints need to be considered: motion constraints and collision constraints.The above will be described in detail next.a)The minimum sum of displacements The minimum sum of displacements is shown in Eq. 13: where there are n segmented teeth T i ði ¼ 1; 2; 3:::; nÞ and m stages for the overall dentition fT n g to move from the initial pose to the target final ideal pose.d ij represents the displacements of the T i at the jth stage and is calculated by Eq. 14: According to the Euclidean distance formula, the displacement of the T i between the j and j À 1 stages can be calculated by the coordinate displacement of the center point, so as to obtain the sum of displacements F 1 of n teeth in m stages.b)The minimum sum of rotations The minimum sum of rotation angles is shown in Eq. 15: where d ij represents the rotation of the ith tooth in the jth stage.As shown in Eq. 16, the rotation of T i between the j and j À 1 stages can be calculated by the vector sum of the rotation angle components, so as to obtain the sum of rotations F 2 of n teeth in m stages.c)Orthodontic movement constraints Only by applying a certain amount of orthodontic force can the teeth keep slow movement by dissolving and rebuilding the alveolar bone (Li et al. 2018).The orthodontic force needs to be within the range of biomechanics.Once the force exceeds the maximum level that the periodontal tissue can bear, the remodeling of the alveolar bone will not be able to keep up with the speed of tooth movement, which will result in orthodontic failure.Therefore, we need to limit the magnitude of the orthodontic force.If all teeth maintain a constant displacement and rotation speed during the orthodontic period, we only need to specify the tooth displacement step d 0 and rotation step d 0 in each stage as shown in Eqs. 17 and 18: where according to stomatological research and clinical experience, d 0 is usually 0.2-0.5 mm, and d 0 is usually 2°* 3°, that is, the displacement value of each tooth in a single stage is at most 0.5 mm, and the maximum rotation value is 3°.d)Orthodontic collision constraints If the path planning scheme is inappropriate, the teeth may still collide and interfere with each other during the movement process.Therefore, we need to carry out collision constraints on the adjacent two teeth at each stage to ensure that there is no collision during the entire movement of the teeth.Let F c be the number of collisions of T i in the whole dentition in a single stage of orthodontics, then F c should satisfy as shown in Eq. 19: where c i is used to judge the collision between T i and its neighbor.As shown in Eq. 20, if there is no collision, c i is 0, otherwise it is 1.
In order to reflect the importance of the collision constraint F c , it is necessary to further add a penalty factor K, which is converted into a new penalty function F 3 as the criterion for collision as shown in Eq. 21: In the formula, the penalty factor K needs to be set larger, the purpose of which is to make F 3 change significantly when adjacent teeth collide.And F 3 is 0 when adjacent teeth do not collide.
To sum up, according to the objectives and constraints in the orthodontic path planning problem, the fitness function SF i of the improved gray wolf optimization algorithm in this paper can be constructed by means of linear weighting as shown in Eq. 22: where i ¼ 1; 2; 3:::; N. N is the size of gray wolf population.k j ðj ¼ 1; 2; 3; 4; 5Þ is the weight coefficient of each item.The choice of k j ðj ¼ 1; 2; 3; 4; 5Þ needs to be handled specifically according to the actual situation of orthodontic path planning.The initial condition of teeth is different from patient to patient, so the orthodontic path planning is also different for each patient.Therefore, the weights of the objective function vary from person to person.The appropriate weights are selected for the corresponding objective functions according to the planning requirements.
It can be known that the smaller the value of SF i , the better the position of the gray wolf.And the gray wolf with the minimum fitness value represents the optimal orthodontic path.

Update the position
The final sub-step is to update the positions of the gray wolves.We use the dynamic weighted position update strategy proposed above to update the positions of the gray wolves in each round of iterations to search for the global optimal solution.The IGWO is shown in algorithm 1.
Several experiments were taken to evaluate the OPP-IGWO method, which was implemented in MATLAB_R2021b and IREACH invisible orthodontic intelligent design system-V2.1.6under MacOS Monterey system, Intel 3.1 GHz dual-core Core i5 and 16 GB memory.

Experimental design and evaluation metrics
This experiment validates our proposed OPP-IGWO orthodontic method.Our dataset consists of tooth models of 30 patients with males (47%) and females (53%) ranging in age from 18 to 40 years old.For each patient, there are two models scanned before and after the treatment.According to Angle's classification, all the three types of malocclusions (Classes I, II, III) are observed in our dataset.
In order to verify that the invisible orthodontic treatment plan generated by our method meets the clinical requirements, the generated orthodontic path needs to be evaluated, which is mainly composed of the following four evaluation indicators: 1. Evaluation index E 1 : corresponding to the objective function F 1 .It is the total displacement of teeth before and after orthodontic treatment (mm).
2. Evaluation index E 2 : corresponding to the objective function F 2 .It is the total amount of rotation before and after orthodontic treatment ( ). 3. Evaluation index E 3 : corresponding to the penalty function F 3 .It is the number of collisions in the whole process of orthodontics.4. Evaluation index E 4 : corresponding to the stage M. It is the number of orthodontic stages required for the whole process of orthodontics.
Weighting these four indicators according to a certain weight, this experiment defines the comprehensive index E score of the orthodontic path effect as shown in Eq. 23:

Comparison with state-of-the-art orthodontic methods
In order to verify the effectiveness of our proposed IGWO algorithm in solving the orthodontic path planning problem, the IGWO proposed in this paper was compared with IGA (Li et al. 2020a, b), Mutil-IPSO (Ma et al. 2021), NSMPSO (Xu et al. 2020), IABC (Li et al. 2020a, b), CGWO (Kohli et al. 2018), HGWODE (Yu et al. 2022), RW-GWO (Gupta et al. 2019) and RLGWO (Qu et al. 2020).The population size of the five intelligent algorithms is set to be 30, the maximum number of iterations is all Each orthodontic path planning method was run 100 times, and the average value was taken as the final result.As shown in Fig. 11, the nine intelligent algorithms show different effects in solving the orthodontic path planning problem. Figure 11a shows an overall presentation of the effects of the nine algorithms, from which it can be seen that the overall effect of our proposed IGWO approach is optimal.With the interpolation initialization strategy, IGWO has a fast convergence speed in the early stage of search.The nonlinear convergence factor and the dynamic weighted position update strategy enable the IGWO to maintain a strong local exploitation capability and continue the search after convergence to a certain extent, thus avoiding falling into local optimum.Figure 11b-i   Table 2 shows the effect evaluation of the orthodontic path generated by different orthodontic methods on patient P's mandibular tooth model data under different evaluation indicators F1-F4.
As can be seen from Table 2, NSMPSO performs best on both E 1 and E 2 , and the orthodontic path planning with NSMPSO has the least total number of positions and total number of rotations of the teeth.However, NSMPSO also has the highest number of collisions, which means that the planned orthodontic paths are difficult to achieve in a real clinical situation.IGWO, HGWODE and RW-GWO have the lowest number of collisions, which indicates that the planned orthodontic paths are fully achievable.RW-GWO also has the lowest number of stages, and its overall performance is close to IGWO, although slightly worse than IGWO.According to the comprehensive index E score , it can be seen that the orthodontic path generated by the IGWO method has the best effect.Compared with IGA, Mutil-IPSO, NSMPSO, IABC, CGWO, HGWODE, RW-GWO and RLGWO, the performance of OPP-IGWO in this paper is improved by 15.46%, 2.24%, 7.18%, 15.99%, 7.67%, 1.01%, 1.42%, 10.23%, respectively.

Comparison with manual orthodontic method
In order to verify the clinical applicability of the IGWO method, in this experiment, the dental models of the 30 patients were divided into three groups according to the degree of malocclusion of the upper and lower jaws, which were recorded as patient groups I, II, and III.Then, we use the IGWO method and the manual planning method to conduct comparative experiments, respectively.Among them, the E 3 collision index is close to 0 for both methods, so it can be ignored.
From Table 3, it can be seen that the total number of stages produced by the IGWO method is significantly lower than produced by the manual method.And the displacements of all orthodontic paths are also generally lower.Compared with manual method, the IGWO method reduced the total amount of tooth displacement before and after orthodontics by an average of 8.32% on the E 1 index.
On the E 2 index, the sum of orthodontic rotations decreased by an average of 5.21%.On the E 4 index, the average number of orthodontic stages was shortened by 34.95%.It is proved that the IGWO method can efficiently and reasonably plan the orthodontic path suitable for invisible orthodontics without considering the influence of clinical complex factors.

Robustness analysis
To verify whether the IGWO proposed in this paper can maintain its robustness in the face of other optimization problems, we selected six commonly used benchmark test functions to test the performance of IGWO against the optimization algorithms mentioned above.The population size of all algorithms was set to 150 and the maximum number of iterations was set to 200.The test functions selected for this experiment are shown in Table 4, where f1 to f3 are single-peak benchmark test functions, with only one global minima, mainly used to detect the local mining ability of the algorithm, including the accuracy size of the solution, convergence speed, etc. f4 to f6 are multi-peak benchmark test functions,  with multiple global minima, mainly used to detect the global search ability of the algorithm, including whether it converges too early, whether it can jump out of the local optimality, etc.To reduce experimental randomness, this experiment runs the above nine optimization algorithms 100 times and takes the mean and standard deviation as the final output results to test the solution accuracy and stability of IGWO algorithm, respectively.The results for the six groups of test functions are shown in Tables 5 and 6.
As can be seen from Table 5, among all the tested functions from f1 to f6, the IGA algorithm has poor performance in finding the best.IGWO, RW-GWO and RLGWO algorithms all show excellent solving ability.IGWO has the highest accuracy at f1, f4, f6.RLGWO has the highest accuracy at f2. RW-GWO has the highest accuracy at f3, f5.
As can be seen from Table 6, the standard deviation of the IGWO algorithm is also on the small side among all the algorithms in all the test functions from f1 to f6.It fully indicates that the IGWO algorithm not only has a strong merit-seeking ability, but also is very stability.
From the accuracy and stability comparison of the algorithm, it can be concluded that the IGWO algorithm with the addition of nonlinear convergence factor and dynamic weight position update strategy better balances local search exploitation and global search capability and demonstrates great superiority in the search accuracy and robustness.
To further evaluate the performance of the improved algorithm visually, this experiment compares and analyzes the convergence speed of these five swarm intelligence optimization algorithms, and Fig. 12 shows the convergence curves and function images of the six groups of test functions f1 to f6.For the convenience of observation, the MATLAB software performs a log transformation of the fitness function to facilitate the differentiation of the convergence curves of different algorithms.
As can be seen from Fig. 12, IGWO performs well in terms of both the merit-seeking ability and the convergence speed.Especially in the multi-peaked functions f4 to f6, the convergence speed and convergence accuracy of IGWO algorithm with the introduction of nonlinear convergence factor and dynamic weighting strategy are significantly higher than other algorithms, especially in the functions f1, f5 and f6.Unlike in orthorectified path planning, the RLGWO has a significant advantage in the merit-seeking ability and convergence speed.It is possible that the reinforcement learning strategy is not applicable to orthotropic path planning.
In summary, the IGWO proposed in this paper performs better in terms of solution accuracy, stability and convergence speed for both single-peaked and multi-peaked functions and satisfies stable convergence with relatively few iterations to complete the search for the global optimal solution.

Conclusion
Aiming at the problem that the current virtual orthodontic system cannot automatically generate invisible orthodontic schemes, this paper proposes an orthodontic path planning method based on improved gray wolf optimization algorithm (OPP-IGWO).Firstly, the tooth model is preprocessed to get the initial tooth model with OBB after segmentation.Then, the target pose of the tooth is determined, and the ideal arch curve is fitted by the Beta curve and Spee curve to obtain the offset of a single tooth in the six degrees of freedom.Finally, according to the initial pose and the target pose of the teeth, the IGWO algorithm is used for orthodontic path planning.The IGWO method is improved mainly from three aspects.(1) Initialization of the gray wolf population based on the initial solution generated by dental interpolation improves the mean value of the initial fitness value of the population and avoids inefficient search of the gray wolf in the initial stage.(2) We introduce a nonlinear convergence factor.The convergence efficiency is slower at the beginning of the iteration, which is beneficial to the global search and exploration of the global optimal solution; at the end of the iteration, the convergence efficiency is accelerated, which is beneficial to the local search and enhances the local exploitation ability.(3) We propose a position update strategy based on a dynamic weighting approach, which introduces a learning rate for each gray wolf to improve the global search ability.And it avoids the situation that the search space is locked in a narrow range because the positions of wolves a, b and d are too extreme.Several groups of experiments show that our OPP-IGWO method is effective in many evaluation indexes, and the planned orthodontic path meets the clinical orthodontic standards.
Although our method can effectively achieve orthodontic path planning, there are still some limitations.1.In this paper, we mainly selected the most representative oval arch morphology and did not consider the special canine round and square arch morphology.And the fitting function used is the more general b-function fitting at present, without considering whether other function fitting is more effective.Therefore, we intend to combine computer vision technology in our future research to automatically identify and analyze the patient's arch morphology and select the appropriate fitting function to fit the ideal arch curve according to different types.2. The orthodontic path planning in this paper does not take into account special situations such as patients' missing teeth and tooth replacement period.Also, by default, the same displacement step and rotation step are used to move simultaneously throughout the treatment, while in the actual orthodontic scenario, the dentist often uses different orthodontic strategies for different teeth at different stages.Therefore, the tooth priority order and the orthodontic force magnitude need to be further investigated.

Fig. 1
Fig. 1 Procedure of invisible orthodontic plan generation

( a )
Feature points of incisors (b) Feature points of fangs (c) Feature points of premolars (d) Features of molars

Fig. 4
Fig. 4 Feature points of different types of teeth

Fig. 5
Fig. 5 Curve fitting method for ideal dental arch

Fig. 6
Fig. 6 Three types of tooth rotation angle

( a )
Before and after teeth arrangement of the upper jaw (b) Before and after teeth arrangement of the lower jaw (c) Before and after teeth arrangement of the whole jaw

Fig. 7
Fig. 7 Visual comparison before and after the row of teeth of patient P

( a )
Initial and target poses of incisors (b) Initial and target poses of the fangs (c) Interpolated pose of incisors (d) Interpolated pose of fangs

Fig. 8
Fig. 8 Initial solution of the orthodontic path

Fig. 9
Fig. 9 Comparison of linear convergence factor and nonlinear convergence factor

Fig. 11
Fig. 11 Iterative diagram of orthodontic path planning shows a separate comparison between IGWO and other algorithms.IGA falls into local optimum earlier.Although Mutil-IPSO converges faster than IGWO in the early stage of the search, it also falls into local optimum at the end, and the final optimization is not as good as IGWO.NSMPSO performs better overall and does not fall into local optimum before the search resources are exhausted.IABC converges slowly.CGWO uses a linear convergence factor and chaotic initialization, converges faster in the early stages of search, but eventually falls into a local optimum.HGWODE and RW-GWO perform similarly to IGWO overall, with slightly worse optimization than IGWO when search resources are exhausted.RLGWO converges more slowly.Overall, IGWO has a faster convergence speed and

Table 2
Evaluation of the orthodontic effect of patient P by different methods

Table 5
Fig. 12 Comparison of the convergence curves of functions f1 to f6