This work addresses the problem of missing tooth crown design using the existing adjacent teeth as shape references, which are shown in Figure 1.

Figure 2 shows the technical process of the proposed algorithm, which mainly consists of the following steps:

(1) Conditional shape model (CSM) construction. Based on a library of clinical oral scan data, the deformable statistical shape models (SSMs) of different teeth were constructed. After that, the inter-teeth shape correlation between adjacent SSMs was modeled using the CSM approach.

SSM (Cootes et al., 1995) parametrically represents the average shape of a class of 3D models and their shape variations, as defined by the following equation, X = X ¯ + ∑ i = 1 p a i Φ i , (1) where, X is the statistical shape model, X ¯ is the average model, Φ i represents the different deformation modes obtained from the principal component analysis (PCA), p is the number of the deformation modes, and a i represents the shape coefficients of the different deformation modes.

To construct the SSM of dental crown, a standard template was aligned to all other clinic oral scan models using the Thin Plate Spline Robust Point Matching (TPS-RPM) (Chui and Rangarajan, 2003) method, to establish the surface vertices correspondence of different templates. In addition, to restrict the vertex distribution, we selected prominent feature landmarks on the occlusal surface of the crown, such as on the grooves, pits, and cusps, based on the anatomical characteristics of different dental positions. These landmarks were used to guide the deformation of the standard template. Figure 3 shows the landmarks selected for different dental positions, where Figures 3A–F are respectively the right maxillary second molar, right maxillary first molar, right maxillary second premolar, right mandibular second molar, right mandibular first molar, and right mandibular second premolar.

From this, we can obtain a set of ideal data, each mesh has the same number of vertices, and the corresponding vertices in different meshes are located in the same anatomical position. N vertices of each mesh can form a vector x , x = x 1 , y 1 , z 1 , . . . , x n , y n , z n T , (2) Where x j , y j , z j is the coordinate of the j -th vertex in the mesh.

To cancel the inter-subject variations in translation, rotation, and scaling, we used Generalized Procrustes analysis (GPA) (Bookstein, 1996) for the sample vectors. Next, we statistically model the aligned shape vectors by principal component analysis (PCA). For the shape vectors represented by Eq. (2), the average model can be obtained by simply averaging the samples as Eq. (3). X ¯ = 1 k ∑ m = 1 k x m , (3) where k is the number of mesh used to set up the SSM.

The singular value decomposition (SVD) of Eq. (4) can be easily calculated to obtain the eigenvectors Φ i and the corresponding eigenvalues λ i , namely, the different modes of deformation and their corresponding variances, with the largest eigenvalue corresponding to the most significant mode of variation. L = x 1 − X ¯ . . . x k − X ¯ , (4) where k is the number of mesh used to set up the SSM.

The parameter a i in Eq. (1) is used to control the shape deformation, which is usually controlled within the range of − 3 λ i , 3 λ i to ensure that the deformation is reasonable.

According to the proposed method, the SSMs are constructed for the missing tooth (described as “son model” in the following) and its adjacent teeth as a whole (described as “father model” in the following). Since the teeth perform important masticatory functions, the size and position of the missing crown mainly depend on its adjacent teeth. In other words, there is a high correlation between the shape coefficients of the two. We parametrically describe this relationship by building the conditional shape model (CSM) (Iglesias and de Bruijne, 2007).

The coefficients of the SSM of the son model S depend heavily on the father model F . The distribution of S can be modeled as a conditional Gaussian distribution P S F under the condition that F is known, P S F = N μ , Σ , (5)

Representations of μ and Σ can be derived according to (Iglesias and de Bruijne, 2007): μ = μ S + Σ S F Σ F − 1 F − μ F , (6) Σ = Σ S − Σ S F Σ F − 1 Σ F S , (7) where, μ is the conditional mean, Σ is the conditional covariance, μ S and μ F are the average shape coefficients of S and F respectively, Σ S and Σ F are the variances of S and F respectively, and Σ S F and Σ F S are the covariance matrices of both.

When a set of SSM distribution coefficients of the father model F is given, the SSM distribution coefficients of the son model S can be obtained according to the above probability model, namely, the son model matching the shape of the father model is obtained.

After constructing the deformable models, in practical application, to obtain the set of shape coefficients of the father model, it is necessary to match the SSM of the father model to the patient’s oral scan model.

According to Eq. (1), the SSM of the father model can be expressed as follows, F = F ¯ + ∑ i = 1 p a i F Φ i F , (8)

The average model vertices of the father model can form a vector F ¯ = x 1 F , y 1 F , z 1 F , . . . , x n F , y n F , z n F T , where x j F , y j F , z j F is the coordinate of the j -th vertex in the mesh.

To improve the computational efficiency, the Iterative Closest Point (ICP) algorithm (Besl and McKay, 1992) was used to initially align the father model with the patient’s oral scan model. The nearest points of F ¯ in the patient oral scan model were searched during each iteration and constituted the vector P = x 1 P , y 1 P , z 1 P , . . . , x n P , y n P , z n P T .

The shape coefficients of the father model a F = a 1 F , a 2 F , . . . , a k F T are computed by solving the linear system to minimize the least squares error, min P − F ¯ − ∑ i = 1 p a i F Φ i F 2 2 , (9)

According to Eq. (6), the shape coefficient of the son model can be expressed as follows, a S = μ S + Σ S F Σ F − 1 a F − μ F , (10)

Then the results of the father model aligned to the patient’s oral scan model can be obtained as follows, F ′ = F ¯ + ∑ i = 1 p a i F Φ i F , (11)

The son model matching with the shape of the father model, S ′ = S ¯ + ∑ i = 1 p a i S Φ i S , (12)

To further improve the alignment accuracy, the TPS-RPM alignment algorithm was used to align F ′ with the patient’s oral scan model again. Since the ICP alignment had matched the two coarsely, TPS-RPM converged at a fast rate to obtain the transformation matrix, which was applied to S ′ .

Finally, the intersection-removing algorithm was used to ensure that there was no unreasonable intersection between the crown and its adjacent teeth. Specifically, the vertices of the crown were looped through, and those located inside the adjacent teeth were moved outward along the direction of the normal vector until all vertices did not need to be adjusted. Finally, the face intersection was detected to ensure reasonable contact between the crown and its adjacent teeth.

The above steps ensure a clear occlusal surface with anatomical features and a harmonious external surface with the morphology of its adjacent teeth (Figure 4A). However, for practical clinical usage, the bottom of the crown must be closely fitted to the preparation.

As shown in Figure 4B, we computed the local curvature for the surface vertices on the preparation model. The vertices with maximum local curvature on the upper part of the model were selected as the candidate points of the finish line (Figure 4C). Next, a polynomial surface is fitted to the candidate points and the crown is cut by this surface to remove the extra parts below the finish line (Figure 4D).

To fit the bottom of the crown to the finish line of preparation, the sign distance function (SDF) was used to control the crown deformation. The SDF is defined as the distance from a point in space to the surface with a positive interior and negative exterior, and the surface is expressed as the 0-equivalent surface of the SDF. First, the mesh of the crown and the preparation were converted into 3D voxels separately and combined into one whole. The SDF values at the vertices of the voxels were calculated, as shown in Figure 4F. In order to create a smooth transition surface from the bottom of the crown to the preparation, the calculated SDF array was diffusion filtered. The cross-sectional view of the crown and the preparation is shown in Figure 4E. To ensure that the shape of the occlusal surface and the preparation below the finish line were maintained, the SDF values in the fixed zone were fixed during the diffusion, while the SDF values in the deformation zone were gradually smoothly and contracted to the finish line. The SDF values after diffusion were shown in Figure 4G. Finally, the 0-equivalent surface of the SDF was converted into a triangular mesh using the marching cube algorithm, resulting in a closed crown that matched the shape of the preparation, as shown in Figure 4H.

Nowadays, additive manufacturing (AM) is increasingly used in the CAM of the dental crown. Compared to the traditional computer numerical control (CNC) technique, AM is more flexible for producing complex structures thus facilitating the manufacturing of complicated support structures inside the object. In this work, skeletonized internal supporting structures of the crown are designed to improve the biomechanical properties and to save the materials for additive manufacturing.

The supporting structure is designed using the triply periodic minimal surfaces (TPMS) method that we previously proposed for generating porous structures (Wang et al., 2022). To effectively represent the porous structures, the period function and the wall-thickness function, which control parameters continuously, were created to control the number of holes (topology change) and the thickness of the structure wall (geometry change), respectively. These two continuous functions, which directly compute integrals and gradients, converted the mechanical problem with given boundary conditions into a computable mathematical model. Moreover, the proposed optimization can be implemented directly on the function representation, only one time of voxel partition is needed to generate the finite elements for integral calculations, instead of multiple times of remeshing of the iteratively updated model. The main optimization process consisted of two steps, where the period function was optimized as a coarse adjustment and the wall-thickness function was used for fine optimization. This two-step optimization can also be accelerated using a global-local radial basis interpolation strategy, which efficiently minimized the compliance of the structures. Finally, a porous structure was obtained with uniform stress distribution.

Specifically, firstly, the initial TPMS-based structure is obtained according to the input crown model and the loading conditions. The initial porosity is 60%, and the wall thickness is set within [0.2, 1.0] mm. Note that the parameters can be adjusted according to the applications. Secondly, a function representation-based optimization formulation is defined to minimize the compliance of the porous shell structure. The bottom of the crown is rigidly fixed to prevent any lateral movement, while disregarding the frictional forces. Additionally, we impose the constraint of vertical pressure exclusively on the occlusal surface of the crown. Thirdly, the aforementioned two-step optimization strategy is exploited to optimize the topology and the geometry, respectively. Finally, we obtain the optimized porous structure of the crown with continuous geometry changes and smooth topology changes that have reasonable stress distribution. In the final processing, an outer shell is added to wrap the porous structures. The default thickness of the shell is set to 0.6–1 mm according to the precision of 3D printers.

Previous to the automatic crown design, the SSMs of the target crown (for the first molars) and its adjacent teeth are constructed. Figure 5 shows the average shape models of the SSMs, where (A-F) are respectively the maxillary second molar, maxillary first molar, maxillary second premolar, mandibular second molar, mandibular first molar, and mandibular second premolar. Figure 5 only shows the models of the right side teeth, whereas the left side teeth models are similar.

It can be seen from Figure 5 that the average shape models maintained the anatomical details of the occlusal surface. Specifically, the maxillary second premolar buccal cusp is similar in height to the lingual cusp, while the mandibular second premolar appears round. The maxillary first molar and the second molar are rhomboidal and generally have four cusps, of which the mesialingual cusp has the largest area as the main functional cusp, but the latter presents a more pronounced rhomboidal shape and a slightly narrower crown. The mandibular first molars typically have five cusps, while the mandibular second molars are cross-shaped. All these detailed features are presented in our shape model.

To observe the shape variation modeling ability of the SSM, different shape coefficients were adjusted to produce deformed shape instances. Figure 6 shows the deformation effects of the first three shape modes. Mode 1 demonstrates the adjustment of width and height aspect ratios. Mode 2 shows the height variation of the lingual and buccal sides of the crown bottom, which is to adapt to the gingival curvature of different patients. Mode 3 controls the relative height changes of the mesial and distal crown and the crest of the tooth.

To construct the CSM, we take the mandibular first molar SSM as the son model, and its five neighboring teeth as a whole SSM as the father model. Then the dependence between the shape coefficients of the two SSMs was constructed as a CSM.

To verify the constructed CSM, we actively changed the shape coefficients of the father model and generated the corresponding son model using CSM. Figure 7 shows the son model shape adaptation following the father model change under the first two modes. Mode 1 mainly focuses on the adjustment of the occlusal position of the father model, and mode 2 mainly shows the change in the tooth height. It can be seen that the son model follows well with the deforming father model in both position and size variations.

Similarly, the CSM of the maxillary first molar was established.

The proposed crown design method was used to test 20 oral scan data randomly selected from our database of maxillary first molar and mandibular first molar, respectively. All the test data were not included in the model training dataset. The first molar and its five adjacent teeth are manually cropped from the organ scan for algorithm validation.

Figure 8A shows the crown design of one representative test data, with an inner surface that matches the shape of the preparation as shown in Figure 8B.

To quantify the accuracy of crown design, we used the manually designed crown by experienced designers as the reference standard. The Chamfer Distance (CD) between our result model and the expert designed to quantitatively assess the accuracy of the SSM matching, C D P 1 , P 2 = 1 P 1 ∑ x ∈ P 1 min y ∈ P 2 x − y 2 + 1 P 2 ∑ y ∈ P 2 min x ∈ P 1 y − x 2 (13) where, P 1 is the point set consisting of the vertices of the crown design, P 2 is the point set consisting of the vertices of the manually designed crown, x and y are the coordinates of the vertices in P 1 and P 2 respectively.

A randomly selected part of the test results is shown in Figure 9. Each column from left to right shows the test data, the deformation results after matching the father model to the test data, the automatic design results of the crown, and the visualization of crown design quality based on CD. For visual assessment, the crown designed by our method has consistent size with the adjacent teeth and complements the surrounding dentition without noticeable gaps or overlaps.

The quantitative test results are shown in Table 1. The American Board of Orthodontics (ABO) accepts crown designs with errors within 0.5 mm (Casko et al., 1998) as a criterion for quality evaluation. In our clinical data experiment, we observed CD of 0.26 ± 0.08 mm for the maxillary first molar and 0.28 ± 0.07 mm for the mandibular first molar.

The internal supporting structure of the crown was optimized using the proposed method. The obtained internal porous structure is shown in Figures 10A, B. The optimized porosity is around 70%, and the wall thickness is adjusted within [0.2, 0.6] mm.

To verify the biomechanical property of the optimized internal supporting structure, finite element simulations were conducted for the designed crown with porous internal structures. The crown material was assumed to be nickel-chromium alloy, with Young’s modulus 214 GPa and Poisson’s ratio 0.3. A nodal load of 0.5 N was applied to the occlusal surface with fixed support conditions to the bottom surface of the model as shown in Figure 11A. Figure 11B shows the distribution maps of the deformation displacement and stress for a representative test data under the given boundary conditions. It can be seen that the stress is uniformly distributed at all internal locations of the crown and the displacement is more concentrated in the cusp region. This result is compared to the simulation of the solid model under the same boundary conditions. The resultant distribution maps of displacement and stress are demonstrated in Figure 11C. The optimized porous internal structure diffused the external pressure and resulted in a more uniform stress distribution than the solid model, and the material used is also saved due to the porous architecture.