BACKGROUND

The objective of an orthodontic treatment planning process is to create a set of valid states that teeth can move between in order to achieve an orthodontically correct final setup. Those states can then be used to design clear tray aligners or other orthodontic appliances. Constraints-based motion between states is a necessary first step in the planning process, as it creates an approximate path that teeth with take between an initial and final setup. It may also be used by a technician in an interactive software system to design a treatment plan or used in a doctor’s office to show the patient an approximation of tooth movement during treatment. Subsequent refinement, however, would likely be required in order to remove collisions between teeth in the interpolated path. A need exists to automatically remove collisions among the states and create a smooth movement path between the states for orthodontic treatment planning or other purposes.

SUMMARY

A first method for automatically removing collisions between digital mesh objects includes receiving a set of digital mesh objects in a first state with one or more collisions between the digital mesh objects, removing the collisions using a perturbation method, and

outputting a set of transformations describing tooth movement between the first state and a second state or tooth positions in the second state.

A second method for automatically removing collisions between digital mesh objects includes receiving a set of digital mesh objects in a first state with one or more collisions between the digital mesh objects, removing the collisions using a mesh deformation method, and outputting the set of digital mesh objects in a second state without collisions between them.

A method for moving digital mesh objects between spatial arrangements includes receiving a set of digital mesh objects in an initial state and receiving a final state for the set of digital mesh objects. The method also includes identifying motion constraints among the digital mesh objects, estimating a number of interpolated states (n) of the digital mesh objects between the initial and final states based upon the constraints, and creating n number of interpolated states for movement of the set of digital mesh objects between the initial and final states.

The states of the digital mesh objects can be used as digital setups for dental or orthodontic treatment planning. BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings are incorporated in and constitute a part of this specification and, together with the description, explain the advantages and principles of the invention. In the drawings,

FIG. 1 is a diagram of a system 10 for removing collisions between digital mesh objects and smoothly move digital mesh objects between spatial arrangements for dental or orthodontic treatment planning;

FIG. 2 is a flow chart of a method to automatically remove collisions between mesh objects;

FIG. 3 is a diagram of a data model used by a genetic algorithm;

FIG. 4 is a diagram illustrating an example of an initial setup and corresponding intermediate and final setups;

FIG. 5 is a flow chart of a method for iterative collision removal of mesh objects;

FIG. 6 is a diagram illustrating mesh objects in collision having multiple collision locations;

FIG. 7 is a flow chart of a method to achieve movement of teeth using a set of constraints on tooth movement; and

FIG. 8 is a diagram illustrating that interpolation between states for single set of constraints can be modified to incorporate different sets of constraints.

DETAILED DESCRIPTION

Embodiments include automated methods to remove collisions between digital mesh objects and smoothly move digital mesh objects between different spatial arrangements. The methods can be used, for example, to generate digital setups for use in dental or orthodontic treatment planning. The methods can also be used with other collections of digital mesh objects in other applications of them.

Methods to automatically generate digital setups for orthodontic appliances are described in US Provisional Patent Application, Serial No. 62/569081, entitled“Automated Process for Intermediate Orthodontic Digital Setup Generation,” and filed October 6, 2017, which is incorporated herein by reference as if fully set forth.

FIG. 1 is a diagram of a system 10 for removing collisions between digital mesh objects and smoothly move digital mesh objects between spatial arrangements

for dental or orthodontic treatment planning. System 10 includes a processor 20 receiving digital 3D models of teeth (12) from intra-oral 3D scans or scans of impressions of teeth, or in other embodiments the system receives user manual input. System 10 can also include an electronic display device 16, such as a liquid crystal display (LCD) device, and an input device 18 for receiving user commands or other information. Systems to generate digital 3D images or models based upon image sets from multiple views are disclosed in U.S. Patent Nos. 7,956,862 and 7,605,817, both of which are incorporated herein by reference as if fully set forth. These systems can use an intra-oral scanner to obtain digital images from multiple views of teeth or other intra oral structures, and those digital images are processed to generate a digital 3D model representing the scanned teeth and gingiva. System 10 can be implemented with, for example, a desktop, notebook, or tablet computer. System 10 can receive the 3D scans locally or remotely via a network.

The 3D scans addressed herein are represented as triangular meshes. The triangular mesh is a common representation of 3D surfaces and has two components. The first component, referred to as the vertices of the mesh, are simply the coordinates of the 3D points that have been reconstructed on the surface - i.e., a point cloud. The second component, the mesh faces, encodes the connections between points on the object and is an efficient way of interpolating between the discrete sample points on the continuous surface. Each face is a triangle defined by three vertices, resulting in a surface that can be represented as a set of small triangular planar patches.

The methods and algorithms described herein can be implemented, for example, in software or firmware for execution by a processor such as processor 20.

AUTOMATICALLY REMOVING COLLISIONS BETWEEN MESH OBJECTS

Described herein are methods that can be used to remove collisions between a set of digital mesh objects, either by perturbing the objects in directions that are likely to remove collisions, or by deforming the shape of the mesh objects. The term digital mesh objects includes mesh objects. A list of the discussed methods and brief description of each method are presented in Tables 1 and

2 below.

Also discussed herein are interactive tools that can be used in an orthodontic treatment planning setting to automatically remove collisions between teeth using these methods.

The description of the methods is divided into two sections. The first section presents a set of intelligent perturbation functions that can be used to shift one or more mesh objects in order to remove collisions without changing their shape. These functions could be applied to remove collisions between a set of teeth without modifying the shape of the tooth mesh. The second section presents techniques that can be used to deform the shape of the mesh objects themselves to remove collisions without shifting their locations. This could be applied to remove collisions between teeth that are falsely in collision due to scanner artifacts. Each approach may be necessary for different parts of the treatment planning process. For example, perturbation may be desirable when attempting to discover a collision-free trajectory between states, whereas mesh deformation may be desirable if collisions are introduced due to artifacts in the mesh that are not present in actual tooth anatomy.

Section 1: Methods to automatically remove collisions using intelligent perturbation functions

I. Directional perturbation

For this approach, an advantageous direction that a mesh or collection of mesh objects are perturbed towards is identified, and the one or more objects are perturbed in this general direction. A flow chart of this approach is shown in FIG. 2.

First, for each object, an advantageous direction, defined as a vector in three-dimensional space, is identified (step 22). This information could either be input by a user or computed based on information such as the vector opposite the penetration between the objects (Method A), the vector perpendicular to the penetration between objects (Method B), or the vector of motion towards an ideal state that does not contain any collisions or contains fewer collisions (Method C). In the next optional step, a vector that is no more than an angle theta from the identified direction is randomly selected for each object (step 24). This optional step provides an element of randomness to the approach that would allow the algorithm to search a space nearby the selected direction for the optimal solution. Next, the magnitude of the perturbation is identified for each object (step 26). The magnitude(s) may be chosen either entirely randomly, through user input, or using some information (e.g., the penetration depth) to define the amount of perturbation, or as a scale factor on the perturbation amount. In the final step, the perturbation is applied to the mesh object(s) by shifting the object(s) by the specified magnitude(s) along their respective vector(s) (step 28).

Method A: Compute the direction of penetration (and penetration depth) between two objects

Aspects of this approach are also referenced in Section 2, Method 1. First, principle component analysis (PCA) is used to fit contact points between two colliding objects to a plane. Other methods for plane fit such as support vector regression can also be used. The penetration direction is identified as the direction normal to the plane (e.g., the vector associated with the smallest principle component), or alternatively, as the mean of the vertex normal directions for all contact points. Penetration depth may also be computed by identifying the vertices that lie within the boundary formed by the contact points (see Section 2, Method 1), then computing the maximum point-to-plane distance when vertices are projected onto the plane along the penetration direction. The sum of the maximum point-to-plane distances for the two colliding objects is the penetration depth.

Method B: Compute the direction perpendicular to the direction of penetration between two objects

The direction perpendicular to the penetration direction may also be computed through PCA or using the plane equation once the plane is computed. Specifically, the vectors associated with the two largest principle components are two directions that are perpendicular to the penetration direction. More generally, the vectors perpendicular to the penetration direction lie on a plane, and any of the vectors may be chosen as the perpendicular direction. Method C: Compute the direction towards an ideal state

In the case that an ideal state that is either collision-free or contains fewer collisions is available, the perturbation direction vector for each object towards this ideal state may be computed as:

II. Propagated direction perturbation

In additional to the objects that are in collision, neighboring objects may also be moved using the perturbation function described above. Perturbing non-colliding objects that are nearby colliding objects has the advantage of creating space for the colliding objects to move into without creating new collisions. The method to perturb colliding mesh objects and their neighbors is similar to the method described above. First, a set of mesh objects that are in the local

neighborhood of the colliding objects is identified. Next, the perturbation direction for each object is identified using one of the methods described above (Method A, B, or C). To minimize unnecessary motion, it may be desirable to set the magnitude of perturbation for each object to a value less than or equal to the magnitude identified for the colliding objects. One option to achieve this outcome is to scale the magnitude of perturbation according to the distance between the colliding objects and the neighboring objects, such that objects that are closer to the colliding objects have larger magnitudes, while objects farther from the colliding objects have smaller magnitudes.

III. Gradient descent

The validity of each point of the trajectory is evaluated using a scoring function, which penalizes teeth collisions. Whenever there is no collision between teeth, this scoring function returns zero, and the score increases with collisions. The gradient descent method finds optimal tooth motion by perturbing their corresponding states in a way that the collision score reduces.

The problem of finding the optimal perturbation that reduces the score function, is a complex nonlinear optimization problem. Although there is no guarantee for finding a global solution, deepest descent methods can be used to find the optimal direction for perturbing states that reduces the collision score to a local minimum.

Since calculating the collision points between the complex geometries such as teeth cannot be represented in a closed mathematical form, the gradient of the score function with respect to object motions should be numerically approximated. To calculate an approximation of gradient, each state can be perturbed many times so that the rate of change of the score with respect to each state can be calculated. The rate of change of the score (i.e., gradient of an objective function) is used to derive a direction vector in which to ultimately move the objects.

First, a subset of the objects in the state are identified by selecting objects that are in collision as well as adjacent objects. For each degree of freedom (or each degree of freedom corresponding to translation) in the degree of freedom vector of each object, the change in the objective function when adding and subtracting a small delta to that degree of freedom is recorded. The resulting vector of changes is normalized and treated as an estimate of the direction of the gradient of the objective function. To create a random perturbation of a state, the negative of the gradient multiplied by a scaling factor is added to the degrees of freedom vectors of the state, along with a small random component. The scaling factor is a random number between 0 and a maximum estimated using the penetration depths of the state.

IV. Simultaneous perturbation stochastic approximation (SPSA)

An alternative method for approximating gradient is called Simultaneous Perturbation Stochastic Approximation (SPSA). In order to approximate the deepest descent direction using SPSA, one can choose a random direction to perturb the current state and calculate the rate of change of the score along that random direction. The main benefit of SPSA method is that in order to have an approximation of the perturbation direction, only two function calls of the score function are enough. In comparison, the numerical gradient approximation requires 2n functions call, where n is the dimension of the state vectors. For the cases that calling the score function is computationally costly, SPSA may result a faster convergence time, since a rough approximation of the deepest descent direction can rapidly progress the optimization toward better regions.

Augmenting SPSA with a momentum technique can also keep track of the history of the estimated directions for a more accurate deepest descent result. In the following section, SPSA method is described in more detail.

Approximating the Steepest Descent Direction using SPSA:

If one considers /{) as a scoring function, and x = [¾, ¾ ...„¾] as the current state, the random direction & = D _¾ ] and the fixed perturbation size c can be used to calculate SPSA of the gradient of / ( ). that is denoted by (:¾£):

In order to avoid numerical instability of calculating each element of D is sampled from a uniform Bernoulli distribution Therefore, the perturbation direction vector

D is a unit norm. In the above definition, the magnitude of the perturbation can be controlled by the perturbation size variable c.

Steepest Descent Optimization Using SPSA:

Let the initial state be denoted by x(Q). the approximation of the steepest descent using the SPSA method would be Therefore, the intelligent perturbation to reduce the score can be calculated by:

where m is the stepsize of the steepest descent. A sequential steepest descent updating method can be used to find a local minima of the score function. In general, at any step the updated state

where ¾ is the stepsize at the i-th iteration. Decreasing the stepsize , during the sequential steepest descent results in a more stable solution, for example /¾ = m/h. where m is the initial stepsize.

Aggregating the History of SPSA Steepest Descent Using a Momentum Method:

Since at each instance of time the SPSA steepest descent encodes the rate of change of the score function along a random direction, the path toward the local minima might be very noisy. The momentum method is a common method for keeping a history of the more recent steepest descent directions for smoothing the steepest descent direction. The momentum direction is the weighted sum of the previously calculated directions where the older directions are receiving lower weight, comparing to the more recent directions. The following equations represent the smoothed direction at the i-th iteration, §{(), calculated using the momentum method

§ CO = C¾-( ) + w§(i - 1),

where w < 1 is the momentum weight, and in this method its value is 0.9.

V. Grid search perturbation

This method performs perturbation by identifying mesh objects that need to be perturbed first using a penetration evaluator, as described earlier in the directed perturbation method.

Additionally, neighboring objects can also be added to this list of possible objects that can be perturbed. The perturbation itself is done by sequentially changing the state of the set of objects, in order to move the objects towards a higher scoring state. For each object, a spatial grid of possible locations is constructed around the object location with a predetermined spatial resolution (e.g. lmm). As an option, some locations can be ruled out using the penetration evaluator to speed-up the search. Then by moving each object to one of the object’s own grid locations, the collision scores can be evaluated. If there are 3 objects and 9 locations for each object, that provides 729 configurations to evaluate. The best scoring configuration is chosen and then this process can be repeated over several iterations as computationally allowed each time reducing the spatial resolution of the grid by a predetermined fraction, for example 0.5. The spatial grid can be constructed in Cartesian coordinates or (for tooth data) archform coordinates. Table 3 provides exemplary pseudocode for implementing this method.

VI. Behavior mode-based perturbation

Some perturb methods could derive higher-level information about the arrangement of mesh objects to choose modes of behavior. For example, if a perturb function detects that an object is in collision with both of its neighbors, it could attempt to create space by flaring those objects through motions that are perpendicular to the penetration directions. These perturb methods might not produce states that immediately improve the value of an objective function, but such states might be critical to producing a path around problem areas.

VII. Genetic algorithm

Described herein is a genetic algorithm (GA) to resolve the collisions for each state (set of teeth). FIG. 3 is a diagram of a data model used by the GA. The GA starts over for each state 34. A GA evolves a population of N chromosomes 36 (solutions) over many iterations (generations), keeping the chromosomes that exhibit high fitness (designated as parents) and deleting the chromosome with low fitness (inferiors). The parents reproduce themselves, and through mutation and crossover operations produce children that pass their successful attributes to the next generation. Evolution proceeds until a criterion is met. In this case, the method ends evolution of an intermediate setup when the method discovers a state which resolves all collisions in that intermediate setup.

Table 4 includes exemplary pseudocode for implementing the GA with the following definitions. A chromosome is a vector of tuples 38, one for each tooth 40, where each tuple 38 consists of a position modifier 42 and an orientation modifier 44.

_ _ Simple Example

FIG. 4 shows a simple example, illustrating the initial setup 46 (malocclusion), an example of one intermediate setup 48, and the final setup 50 for a simplified arch of 4 teeth. The GA is run on this intermediate state in order to resolve the collision (see the middle two ellipses in setup 48). In this case, the algorithm might create a population of N = 100 chromosomes, where each chromosome comprises a vector of tuples, where each tuple comprises a small displacement and a small change in orientation (e.g., <0.04 mm, 0.09mm, 0.01mm>; 1.5 degrees along a given axis). The algorithm iterates the following for as many generations as necessary in order to produce a chromosome without collisions: A) measure the fitness of each chromosome (e.g., as a count of collisions), B) shuffle the population into groups of four and flag to the two lowest-fitness chromosomes for deletion, and C) replace the inferior chromosomes with copies of the parents that have undergone crossover and mutation. Iterate until done. Eventually, after many generations, the teeth in the intermediate state have been moved such that the teeth no longer have collisions. Then the process is restarted, using the next intermediate setup as input.

Interactive Tools for Orthodontic Treatment Planning from Perturb Methods

In addition to applications in algorithms that automatically remove collisions in a series of mouth states, methods that perturb teeth toward a more valid state have many potential applications in interactive tools for orthodontic treatment planning, including interactively generating intermediate and final setups, and for dental treatment planning, including planning optimal placement of an implant.

In preparing a specific arrangement of teeth for orthodontic treatment, a technician may be trying to achieve certain qualities, such as aesthetic appearance or directives from an orthodontist. However, it may be difficult to manually manipulate teeth at a small scale, so these goals are satisfied without creating collisions between the teeth. Even if the treatment planning software does not allow the teeth to be moved into collision, extra work may be required to manipulate teeth (either one at a time or in groups) into the desired state without creating collisions.

An interactive tool that allowed a perturb function to be applied to reduce or eliminate small collisions would allow the technician to focus less on collisions and more about other goals. The technician could place teeth into a state that satisfies their high-level goals but has collisions, and use the tool to find an extremely similar state that does not have collisions.

Many perturb functions described above rely on scale factors that place bounds on the distance from the original state to the perturbed state. These scale factors could be linked to sliders or other controls in a graphical user interface that allow the technician to control how extreme the motions of the teeth are as collisions are resolved.

Embodiments of this invention could also be applied in a dentistry setting to enable automated or semi-automated planning of the location of a crown, bridge, or implant. In this application, a dentist could roughly place an implant in a digital arch form such that it achieves general aesthetic goals. Intelligent perturbation functions could then be used to automatically minimize or removing both collisions between teeth on a single arch and collisions between teeth on opposite arches. Mesh deformations (described in the next section) could also be applied to the implant in order to design an implant that fits correctly in a patient’s mouth. As an alternative to collisions, another score such as penetration depth or number of colliding points could also be used.

Section 2: Methods to automatically remove collisions through mesh deformation

Various parts of the digital orthodontic workflow may produce valid arrangements of teeth that contain collisions. There are two main causes for these collisions. First, while the teeth in their maloccluded state will not contain collisions, related geometry (such as roots) with collisions may be automatically generated. Second, technicians may intentionally create arrangements of teeth with small collisions between teeth or roots. Their rationale is that these collisions are between portions of the mesh that may represent scanning issues or segmentation errors, so they should be disregarded.

Since the various processes that optimize arrangements of teeth assess validity by checking for collisions between teeth, such collisions (especially in the final setup) are problematic. However, automatic methods may be used to eliminate collisions by modifying the geometry. These methods could also be applied interactively so that technician could confirm that the modifications to the geometry are not too extreme or are only applied to regions or objects that the technician selects.

Method 1: Automatically remove collisions between meshes through local mesh deformation

Described herein is an iterative approach for collision removal whereby the portions of mesh objects that are in collision are identified and shifted until collisions are removed.

The process for iterative collision removal is shown in the flow chart in FIG. 5. First, the contact points between two colliding objects are identified (step 52). It is possible that objects in collision may have multiple collision locations, as illustrated in FIG. 6 showing collision locations 64. Therefore, the next step uses minimum linkage clustering to identify the contact points associated with each collision location (step 54). For each collision location, the contact points are fitted to a plane using PCA (step 56). The penetration direction may be defined either as the vector that is normal to the fitted plane (i.e., the vector associated with the smallest principle component) or the mean of the vertex normal directions at the contact points (step 58). Two ways are possible to identify mesh vertices that are in collision (step 60). In the first approach, mesh objects are cut at the plane identified by PCA, and only the sub-meshes that are closest to one another are retained. The vertices of both sub-meshes are then projected onto the plane, and vertices that are within the convex hull of the contact points are identified as being in collision. In the second approach, vertices in opposing meshes that are within a small distance, d, of each other when projected onto the plane are identified as being in collision. Colliding vertices on each mesh object are subsequently shifted by a small amount, d, in a direction opposite to the penetration direction (step 62). This process is repeated iteratively until no collisions remain between the two objects.

Method 2: Global mesh shrinking

Table 5 provides exemplary pseudocode for implementing a global mesh shrinking method.

The advantage of this method is that it preserves the overall shape of the meshes being modified and does not introduce sharp edges into previously smooth areas. However, if large shrink amounts are required to remove collisions, a large amount of error may be introduced.

This method may be better suited for automatically generated geometry, such as roots, which should retain smooth surfaces while not colliding with each other or with crown meshes.

SMOOTHLY MOVING MESH OBJECTS BETWEEN SPATIAL ARRANGEMENTS

Described herein is an approach for smoothly moving a set of digital 3D mesh objects between two spatial arrangements. This method incorporates motion constraints on the per-state movement of each object, allowing for multiple objects to be moved at different rates according to their constraints. In an orthodontic treatment planning process, this constraints-based approach could be used to model tooth movement between states according to the amount of force applied to each tooth.

Orthodontic treatment planning requires moving a set of teeth from one stage to the next until arriving at the final setup. A basic approach to achieve this goal is to move teeth between stages using linear interpolation. Basic interpolation, however, assumes that all teeth move at a uniform rate, which is not the case. Rather, the amount of tooth movement for each tooth is based on the amount of force applied to the tooth (determined by the treatment appliance) and the anatomy (tooth shape, root length, jaw thickness, and others). To achieve realistic movement of teeth, one approach uses a set of constraints on tooth movement to reflect the amount of motion that is achievable on a per-tooth basis, as illustrated in the flow chart of FIG. 7.

First, the motion constraints per tooth are identified (step 66). These constraints may be provided manually by a clinician, established through heuristics, or computed using biomechanical modeling. Constraints may also be informed by other modalities such as x-ray or computed tomography (CT), which can be used to determine the root shape, root length, and jaw thickness. This information would then be combined with the treatment appliance information to determine constraints on tooth movement. Once the constraints have been established, the total number of interpolated states is estimated by subdividing total tooth movement by the constraints (steps 68). To subdivide at a higher resolution than the constraints, the constraints may be multiplied by a scale factor. Table 6 provides exemplary pseudocode for implementing these steps 66 and 68.

Tooth movement may be subdivided into the number of states, n, using two approaches. In the first approach (“same finish”), tooth movement is planned such that all teeth finish moving at the same time (step 70). This means that teeth with the same constraints but a different total amount of movement will move at different rates so that they finish moving at the same time. A same finish approach could be used to minimize patient discomfort during treatment. In the second approach (“fast finish”), tooth movement is planned such that each tooth finishes moving as quickly as is permitted by the constraints (step 72). This means that teeth with the same constraints but a different total amount of required movement will move at the same rate and finish moving at different times. A fast finish approach could be used to improve aesthetic appearance more quickly during treatment.

Same finish:

Tooth translation is subdivided into the number of states, n, by dividing the total translation amount evenly among the states. Rotation is subdivided into the number of states, n, using spherical linear interpolation (SLERP) of quaternions. SLERP has several advantages over linear interpolation: The shortest path between two angles is guaranteed; and the velocity of the object remains constant throughout the rotation.

Fast finish:

A tooth-specific number of states, , is computed for each tooth. For each tooth, tooth translation is subdivided into the number of states, , by dividing the total translation amount evenly among states. Rotation is subdivided into the number of states, , using SLERP. The interpolation for each tooth is subsequently zero-padded so that the total number of states is n.

Enabling options for more orthodontist control and/or hybrid treatment

By controlling the constraints (and thus, resolution) in the interpolation parameters and/or staging thresholds, either a coarse or fine orthodontic correction can be accomplished. These constraints can be dictated by the amount of control an orthodontist would like to have in the treatment process. High resolution (small constraints) give an orthodontist greater ability to control and monitor the treatment as it progresses, since it results in a fine orthodontic correction that could imply more visits to the orthodontist office and smaller, more precise motions exerted by the clear tray aligners (CTAs). Low resolution (large constraints) provide less control, which could be desirable for less complex cases and to keep the costs low. This multiplicity of options can also translate to a hybrid treatment process. During the larger movement phase, traditional brackets and wires can be involved to result in a stage at which CTAs can take over.

The method above could be adjusted to allow for forces that change throughout the course of treatment, as illustrated in FIG. 8, allowing for treatment planning with a set of different appliances (e.g., brackets and wires followed by clear tray aligners). To modify the forces, first define an intermediate state where treatment modality will be switched. Next, define the constraints on tooth movement from the initial state to the intermediate state (cl) and from the intermediate state to the final state, (c2), and perform interpolation. FIG. 8 is a diagram illustrating that interpolation between states for single set of constraints 74 can be modified to incorporate different sets of constraints 76 due to, for example, a change in the corresponding treatment appliance.

This process is also applicable to CTAs made from different materials which have different inherent forces but at varying costs. This will allow orthodontists to communicate to the patient about different treatment options based on cost (material + time cost + office visits), treatment time, aesthetics (CTA time), and other factors.