{'fraction_non_alphanumeric': 0.0268899634390848, 'fraction_numerical': 0.009133676665225196, 'mean_word_length': 5.052664974619289, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Open-Full-Jaw: An open-access dataset and pipeline for finite element models of human jaw

Developing computational models of the human jaw acquired from cone-beam computed tomography (CBCT) scans is time-consuming and labor-intensive. Besides, a quantitative comparison is not attainable in the literature due to the involved manual tasks and the lack of surface/volumetric meshes. We share an open-access repository of 17 patient-specific finite-element (FE) models of human jaws acquired from CBCT scans and the utilized pipeline for generating them. The proposed pipeline minimizes model generation time and potential errors caused by human interventions. It gets dense surface meshes and provides reduced conformal surface/volumetric meshes suitable for FE analysis. We have quantified the geometrical variations of developed models and assessed models' accuracy from different aspects; (1) the maximum deviations from the input meshes, (2) the mesh quality, and (3) the simulation results. Our results indicate that the developed computational models are precise and have quality meshes suitable for various FE scenarios. Therefore, we believe this dataset will pave the way for future population studies.


Introduction
Finite element modeling (FEM) is a numerical approach for predicting responses of different tissues under physical loads, which can be difficult or impossible different tissues (e.g., periodontal ligament, gingiva, and alveolar bone) or orthodontic appliances [3]. Except for a few recent studies [4,5,6], almost all of the previous studies in the field are limited to single-patient analysis [7,8,9,2,10,11], in which the results might not be generalized to a larger population with high geometrical variations in the teeth, periodontal ligament (PDL), and bone anatomies [4,6].
The main reason for using a single model in the literature is that developing accurate computational models of the human jaw is challenging and involves timeconsuming and labor-intensive processes such as segmentation, geometry reconstruction, geometry processing, remeshing, and mesh simplification tasks. For instance, generating a complex and highly detailed finite element (FE) model of the entire human jaw can take up to several months per scan [10]. Therefore, developing several patient-specific FE models may not be feasible for many researchers. In addition, currently, there are no publicly available datasets of full dentition human jaw to be used by researchers, except for two studies [10,17] with a limited number of studied subjects.
In other words, in almost all of the studies focusing on full dentition or single/multiple tooth analysis [2,11,7,8,6,9,2,14,15,16], the utilized geometries, volumetric meshes, and FE models have not been made publicly available, which makes it difficult to reproduce and compare the results. Table 1 presents an overview of the related studies in the literature by providing details on the studied cohort, discretization type, and availability of the models.
As one of the few studies with public data, the OpenMandible [10] provides detailed geometries of one mandible structure and all teeth obtained from a dried male skull. The study scans the mandible in two different steps to provide detailed teeth structures, i.e., pulp and enamel. First, the mandible was scanned using a CBCT scanner with a voxel size of 0.133 mm. Second, all mandibular teeth were removed from the bone before being scanned by the micro-computed tomography (micro-CT or µCT). Micro-CTs are high-resolution CT scans that are normally acquired from dead specimens or cadavers due to high x-ray exposures. However, the detailed geometries obtained from a single mandible cannot cover geometrical variations across different patients in the population.
The recently introduced OpenJaw Dataset [17] provides open-access reconstructed geometries of three patients' mandibles acquired from CBCT scans. Each patient's data includes surface meshes of the reconstructed mandible, teeth, and PDL geometries. In the geometry reconstruction step, the utilized scans in the study (one with 0.15 mm and two with 0.3 mm voxel sizes) were upsampled to the same resolution of 0.15 mm [5]. Although this study provides more public samples, generalizability to the population remains a problem. Moreover, both of the abovementioned studies include manual tasks in different geometry-processing or meshing tools making it difficult for other researchers to reproduce their meshes using the unprocessed reconstructed meshes.
Our main contributions in the Open-Full-Jaw study can be summarized as, 1. We provide an open-access dataset of different patient-specific models of the human jaw, including the maxilla, mandible, full dentition, and the PDL geometries obtained from CBCT scans of 17 patients. It is the largest publicly available dataset with validated segmented geometries and quality volumetric meshes that can directly be used in FEM studies. Furthermore, to the best of our knowledge, it is the first repository containing the maxillary jaws of different patients. 2. We introduce a unique repository (https:// github.com/diku-dk/Open-Full-Jaw) containing (1) clinically validated segmented geometries and the resulting dense irregular surface meshes; (2) the quality and adaptive volumetric/surface meshes to be used in FE simulations; (3) the automatically generated FEM files for tipping and biting scenarios used for the FE analysis of this work; (4) the principal axes of every patient's tooth providing great information for the users to automatically set up different loading conditions. 3. For reproducibility, we share our pipeline developed based on open-source meshing tools [18,19] to generate the models of this study. This python-based library automates the FE model generation process, including geometry processing and re-meshing tasks with minimal human intervention, by setting a few required parameters. 4. This pipeline allows other researchers in the field to generate quality volumetric meshes and FE models directly from dense and uncleaned meshes with minimal human intervention. This will help other researchers to easily extend their datasets without spending much time and effort on manually cleaning up the meshes and non-trivially producing conformal meshes. 5. Our pipeline ensures conformal meshes in the contacting interfaces without any undesired gaps or pen-  [10]; CG: Congruent interface (no gap/penetration in the contacting interfaces); Non-CG: Non-congruent interface; CM: Conformal Mesh; Non-CM: Non-Cnformal Mesh. Obtained from cadaver or dried skull. * See Table 2 for further details., † Non-smooth top surface of PDL, see [12,13,15] for illustrations.
etrations and provides adaptive meshes that are vital for reducing the total number of elements while using finer meshes in specific regions, e.g., teeth sockets, alveolar crest, and alveolar process.
All in all, we believe the Open-Full-Jaw dataset can be used for various intra-and inter-patient analyses [4,6] such as intact tooth movement modeling, bite force estimation, restorative procedures modeling including cavity fillings and dental implants, just to name a few. Besides, it can greatly impact the reproducibility of future studies. For the reproducibility of this study, we use open-source meshing tools. Still, the reconstructed geometries provided in our dataset can be imported into any desired open-source or commercial meshing tools or FE frameworks to re-mesh and generate computational models.

Important traits required for a successful FEM
Developing a patient-specific FE model begins with segmenting/annotating the desired regions of the medical scan obtained from the patients. Next, the segmented regions are reconstructed as surface meshes generally composed of irregular dense meshes with no guarantees of manifoldness, watertightness, or absence of selfintersection, which are crucial for developing stable and accurate computational models. Hence, one needs to generate quality meshes from the exported dense and irregular meshes that are not necessarily guaranteed to have the mentioned criteria.
Moreover, different preprocessing steps such as geometry processing, mesh reduction, and re-meshing are es-sential for developing FE models from image-based reconstructed geometries. When modeling geometries with shared contacting interfaces, each of the mentioned processes can produce errors on the contacting surfaces and result in undesired gaps/penetrations between them. In the cases where two adjacent segments are watertight, it is still challenging to discretize the segmented domains such that they agree on the same discretization on the shared contacting interfaces. Therefore, the focus of this section is on essential aspects needed to be considered for the discretization of the computational domains with shared contacting interfaces; we also discuss potential options for developing a proper FE model of the human jaw.

Congruent contacting interfaces
The discretization of the computational domain is an essential step for developing computational models considering the biomechanical behavior of tissues having sliding or bilateral contact with other domains/tissues. It mainly affects the numerical stability of the computational models when soft structures, e.g., PDL, contact hard tissues like bone or teeth. Besides, a coarse or poor-quality discretization can itself cause a locking effect and influence the accuracy of the stress/strain concentrations compared to that of the values measured in vivo.
When two contacting domains are reconstructed, discretized, or re-meshed independently, undesired gaps/penetrations are inevitable, causing two different boundary definitions for an identical contacting interface between the two adjacent/contacting domains. The agreement of the two contacting geometries on the contacting boundary or shared interface can be analyzed in two different levels; first, on the geometry level and based on the curvature of the contacting surfaces; second, on the Figure 1: A visualized comparison between the ε-congruent and conformal meshes. A: the ε-congruent contacting surfaces with gaps/penetrations in the contacting interfaces. B: the same structure with the multi-domain discretization resulting in conformal meshes. C: a close-up cross-section view of the gaps/penetrations at the root level. D: a schematic illustration of the ε-congruent and conformal meshes.
mesh-based level focusing on the agreement on the identical discretization of the contacting interface in terms of the position of vertices, edges, and faces. The former is called here the "interface congruency". In the biomechanical field of study, different theoretical and computational studies analyze the effect of the "ball-and-socket" joint congruencies such as in shoulder, hip, and temporomandibular joints to analyze their instabilities and dislocations under different circumstances [20,21]. It should be noted that the current study only focuses on developing computational models of human mandible and maxilla for tooth movements and the mesh congruency of the contacting surfaces instead of congruency of the "ball-andsocket" joints, and the utilized "congruency" term needs to be distinguished from those in biomechanical studies for analyzing the joint congruencies [20,21].
As mentioned before, by using the conventional oneby-one domain discretization, mesh reduction, or quality meshing processes, it is challenging to achieve congruent surfaces due to the generation of gaps/penetrations between two contacting surfaces (see Figure 1). We propose using signed distance fields in the contacting surfaces to quantitatively evaluate the error between two contacting regions. The signed distance function, φ Ω (x), for a domain Ω ⊂ R 3 and an arbitrary point x ∈ R 3 is defined as where ∂Ω denotes the boundary of the domain Ω, Ω c represents the complement of Ω, and d(x, ∂Ω) indicates the Euclidean distance between the point x and the boundary of the domain ∂Ω. We measure the congruency of the surfaces by calculating the signed distance field of each contacting surface with respect to the other domain as where ∂Ω 1 and ∂Ω 2 denote the boundaries of the domains Ω 1 and Ω 2 , respectively, and ∂Ω 1c and ∂Ω 2c refer to the contacting surfaces of the domains. The computed ε values for two contacting domains are then called epsiloncongruent surfaces. The epsilon value indicates the measurement error, where the values close to zero refer to completely congruent surfaces.

Conformal mesh interfaces
A special case for the congruent surfaces is "conformal meshes" or "mesh conformity" and can be described as the identical discretization of the contacting interfaces. That is to say, for two contacting domains/geometries like Ω 1 and Ω 2 , the contacting surfaces ∂Ω 1c and ∂Ω 2c are assumed identical. More specifically, they share identical vertices, edges, and elements on the contacting interfaces. In general, applying the conformal mesh criterion on multi-domains with contacting surfaces is a challenging process, as most of the geometry and mesh processing algorithms used in different free software produce the quality meshes per domain independently. Therefore, the contacting surfaces should be identified and combined as a different step to create identical contacting interfaces on the mesh of each domain.
Finally, the volumetric meshes can be generated by enforcing the meshing algorithm to fully preserve the input surface mesh. This raises questions about the final  mesh quality, as the additional surface preservation constraint affects the mesh generation and the optimization algorithms. Hence, the best approach for generating conformal meshes is to consider all contacting domains simultaneously while generating the volumetric meshes. Generating multi-domain volumetric meshes would generate zero-congruent and conformal meshes that omit any numerical errors caused by non-congruent contacting meshes. Therefore, it is essential to utilize a meshing algorithm that generates volumetric meshes considering multidomain boundaries and performing boolean operations on the input surface meshes by using an implicit representation of the input meshes. In this study, we use fTetWild [18] that supports all aspects mentioned above for multidomain volumetric meshes. This provides a volumetric mesh of teeth-PDL-bone geometries in which there are shared points in the teeth-PDL and PDL-bone contacting surfaces.

Proper modeling of the PDL layer
PDL is a thin structure that connects alveolar bone to the teeth cementum and acts as a shock absorber in the chewing or biting process. It also plays an important role in transferring load from teeth to the bone in orthodon-tic treatments; when triggered with enough orthodontic forces, it results in a bone remodeling process. As this thin and soft structure shares interfaces with two hard tissues of the teeth and bone geometries, the computational model may become unstable if the PDL geometry does not have congruent surfaces with teeth and bone in the contacting interfaces.

Quality mesh generation
Tetrahedralization algorithms usually require a set of closed, self-intersection free surface meshes as input, which is a requirement that is challenging to guarantee in our setting without manual interaction. In fact, traditional medical imaging pipeline require the large majority of the manual processing in this cleanup phase, as selfintersection, holes, or other imperfections normally appear when processing real-world data, and are especially common (and time consuming) when dealing with thin layers, like the PDL layer.
We propose a very different approach, where we accept that these imperfection exist, and design a meshing pipeline that tolerates and automatically heals them. We base our approach on the method recently introduced in [22,23]: instead of meshing one domain at a time, TetWild 5 Figure 2: Reconstructed 3D geometries from 17 CBCT scans. Note that the maxillary jaw is removed in some scans due to having several metal artifacts and/or missing teeth. The anatomical variations in different jaws and teeth of our dataset indicate the necessity for introducing such a model repository for assessing the generalizability of the FEM results in the related clinical/population studies. mesh quality, as the additional surface preservation constraint affects the mesh generation and the optimization algorithms. Hence, the best approach for generating conformal meshes is to consider all contacting domains simultaneously while generating the volumetric meshes.
Generating multi-domain volumetric meshes would generate zero-congruent and conformal meshes that omit any numerical errors caused by non-congruent contacting meshes. Therefore, it is essential to utilize a meshing algorithm that generates volumetric meshes considering multi-domain boundaries and performing boolean operations on the input surface meshes by using an implicit representation of the input meshes. In this study, we use fTetWild [18] that supports all aspects mentioned above for multi-domain volumetric meshes. This provides a volumetric mesh of teeth-PDL-bone geometries in which there are shared points in the teeth-PDL and PDL-bone contacting surfaces.

Proper modeling of the PDL layer
PDL is a thin structure that connects alveolar bone to the teeth cementum and acts as a shock absorber in the chewing or biting process. It also plays an important role in transferring load from teeth to the bone in orthodon-tic treatments; when triggered with enough orthodontic forces, it results in a bone remodeling process. As this thin and soft structure shares interfaces with two hard tissues of the teeth and bone geometries, the computational model may become unstable if the PDL geometry does not have congruent surfaces with teeth and bone in the contacting interfaces.

Quality mesh generation
Tetrahedralization algorithms usually require a set of closed, self-intersection free surface meshes as input, which is a requirement that is challenging to guarantee in our setting without manual interaction. In fact, traditional medical imaging pipeline require the large majority of the manual processing in this cleanup phase, as selfintersection, holes, or other imperfections normally appear when processing real-world data, and are especially common (and time consuming) when dealing with thin layers, like the PDL layer.
We propose a very different approach, where we accept that these imperfection exist, and design a meshing pipeline that tolerates and automatically heals them. We base our approach on the method recently introduced in [22,23]: instead of meshing one domain at a time, Table 2: Specifications of the utilized CBCT scans. Different scans with various filed-of-views and slice thicknesses are used in this study. All scans are converted to an identical slice thickness of 0.15 mm to avoid biases in the geometry reconstruction process, especially when applying smoothing filters to eliminate noises on the segmented teeth and bone geometries. Note that the maxillary jaws of the patients with many metal artifacts and/or missing teeth are removed. The scan obtained from 3DSlicer's "Sample Data" module, titled "CBCT-MRI Head". † The mandible includes partially erupted wisdom tooth/teeth. * The maxilla includes an impacted wisdom tooth.
TetWild meshes the entire volume of the bounding box containing the soup of triangles of all surfaces of interest. The triangles are approximated with faces of the tetrahedral mesh. This procedure does not require clean input geometry, it can tolerate degenerate input triangles, selfintersections, and holes [22]. In the original TetWild algorithms, the final mesh is extracted using a robust filtering procedure based on either flood fill or the generalized winding number [24], in both cases assuming that the user is interested in a single material mesh. We propose a novel filtering procedure for the filtering of multi-material tetrahedral meshes common in medical imaging in Section 3.5.5, and we show that our extension of TetWild is ideal for constructing our dataset, as it removes the expensive and tedious manual cleanup of the input surface meshes.

Controlling shared interfaces using volume mesh generation
This section reviews the entire process for developing the FE models of 17 jaws presented in Figure 2. To be more specific, Section 3.1 describes the utilized criteria for the cohort selection process; Sections 3.2 and 3.3 provide a detailed description for the CBCT segmentation and clinical validation; the geometrical variations of the reconstructed mandibles and maxillae are investigated in Section 3.4 based on widely used clinical landmarks; and finally, the different steps of the used pipeline for developing high-quality volumetric meshes of human jaws are studied in Section 3.5.

Cohort selection
We use CBCT scans of 17 different patients with different voxel resolutions from 3Shape A/S in-house dataset. Various criteria including the original voxel size of the scan, minimal metal filling artifacts, and the absence of implants and severe periodontal diseases are considered in selecting the mentioned cohort. Also, the selected cohort has no evidence of maxillofacial surgery or skeletal diseases. Sensitive information of the patients such as name, age, and gender are stripped due to the General Data Protection Regulation (GDPR) rules. The utilized scans were acquired from the patients by their associated doctors/orthodontists as a part of treatment plans, and the authors had no roles in the acquisition process.

Data specifications and geometry reconstruction
To reconstruct the patient-specific geometries, first, the scans are imported in 3DSlicer [30] in the standard Digital Imaging and Communications in Medicine (DICOM) format. Table 2 provides details of the scans utilized in this study. Next, according to the pre-evaluation criteria (metal fillings or implant artifacts), we decide on which jaws are suitable to be segmented from the scan.   [25,26,27,28,29], as illustrated in Figure 3. The perpendicular distance from the left canine eminence (CE l ) to the alveolar crest SCT Angle The angle between the lines connecting the left canine eminence (CE l ) to the anterior nasal spine (S ) and the left maxillary tuberosity (T l ) TST Angle Anterior nasal spine angle: the angle between the lines connecting the anterior nasal spine (S ) to the left and right maxillary tuberosity (T l , T r ) TT Width Posterior maxillary width; the distance between the left and right maxillary tuberosity (T l , T r ), inspired from [27] The resolution of the selected CBCT scans is at most 0.3 mm since based on our experience, accurate tooth-bone segmentation of the scans with slice thicknesses above 0.3 mm is very challenging. Besides, applying smoothing functions with identical kernel sizes results in smoother segments for coarse voxel sizes, making it difficult to remove the segmentation noise (e.g., rugged surfaces) while preserving the fine details in desired regions (e.g., the alveolar crests of teeth sockets and root apexes). As a result, to avoid biases in the geometry reconstruction process, we convert all scans to an identical resolution of 0.15 mm, as summarized in Table 2, by oversampling a cropped region of interest (ROI) containing the jaws using linear interpolation.
The tooth-bone segmentation is performed based on the semi-automatic watershed algorithm [31] provided in 3DSlicer's SegmentEditorExtraEffects extension. The result of the watershed algorithm is then refined to correct the misclassified teeth and bone. Later, the geometries are smoothened using the 3DSlicer's standard median and joint smoothing modules [32,33]. The joint smoothing method [33] smooths the adjacent segments simultaneously and enforces watertight interfaces between them. The tooth-bone segmentation procedure proceeds until the segmentation accuracy meets our clinical validation criteria (see Section 3.3). Finally, the segmented regions are exported as surface meshes in the Object file format (OBJ).

Clinical validation of the teeth-bone segmentation
The segmented teeth-bone geometries of all patients are validated by clinical experts. This is done by identifying the existence of any periodontal diseases and categorizing each patient into one of the no, mild, or moderate periodontal disease categories. The patients with severe periodontal diseases are excluded from the analyses based on the cohort selection criteria mentioned in Section 3.1. Moreover, we assess the accuracy of segmentation in areas close to teeth/bone borders, cervical regions of bone around teeth sockets, tooth-bone interfaces, and roots.

Geometrical variations of the dataset
We evaluate the geometrical variation of the reconstructed geometries using 3D morphological and cephalometric landmarks and measurements adapted or inspired from the literature [25,26,27,28,29], as illustrated in Figure 3 and described in Table 3. Accordingly, the measured values per patient and the mean and standard deviation of each morphological measurement across all patients are listed in Table 4. Note that even though the standard deviation of each measurement may seem rather small, such 7 Figure 3: The utilized morphological landmarks for analyzing and quantification of the geometrical variation of the reconstructed mandibles and maxillae in different models. Table 3: Description of the utilized morphological landmarks obtained based on the literature [25,26,27,28,29], as illustrated in Figure 3. The perpendicular distance from the left canine eminence (CE l ) to the alveolar crest SCT Angle The angle between the lines connecting the left canine eminence (CE l ) to the anterior nasal spine (S ) and the left maxillary tuberosity (T l ) TST Angle Anterior nasal spine angle: the angle between the lines connecting the anterior nasal spine (S ) to the left and right maxillary tuberosity (T l , T r ) TT Width Posterior maxillary width; the distance between the left and right maxillary tuberosity (T l , T r ), inspired from [27] The resolution of the selected CBCT scans is at most 0.3 mm since based on our experience, accurate toothbone segmentation of the scans with slice thicknesses above 0.3 mm is very challenging. Besides, applying smoothing functions with identical kernel sizes results in smoother segments for coarse voxel sizes, making it difficult to remove the segmentation noise (e.g., rugged surfaces) while preserving the fine details in desired regions (e.g., the alveolar crests of teeth sockets and root apexes). As a result, to avoid biases in the geometry reconstruction process, we convert all scans to an identical resolution of 0.15 mm, as summarized in Table 2, by oversampling a cropped region of interest (ROI) containing the jaws using linear interpolation.
The tooth-bone segmentation is performed based on the semi-automatic watershed algorithm [31] provided in 3DSlicer's SegmentEditorExtraEffects extension. The result of the watershed algorithm is then refined to correct the misclassified teeth and bone. Later, the geometries are smoothened using the 3DSlicer's standard median and joint smoothing modules [32,33]. The joint smoothing method [33] smooths the adjacent segments simultaneously and enforces watertight interfaces between them. The tooth-bone segmentation procedure proceeds until the segmentation accuracy meets our clinical validation criteria (see Section 3.3). Finally, the segmented regions are exported as surface meshes in the Object file format (OBJ).

Clinical validation of the teeth-bone segmentation
The segmented teeth-bone geometries of all patients are validated by clinical experts. This is done by identifying the existence of any periodontal diseases and categorizing each patient into one of the no, mild, or moderate periodontal disease categories. The patients with severe periodontal diseases are excluded from the analyses based on the cohort selection criteria mentioned in Section 3.1. Moreover, we assess the accuracy of segmentation in areas close to teeth/bone borders, cervical regions of bone around teeth sockets, tooth-bone interfaces, and roots.

Geometrical variations of the dataset
We evaluate the geometrical variation of the reconstructed geometries using 3D morphological and cephalometric landmarks and measurements adapted or inspired from the literature [25,26,27,28,29], as illustrated in Figure 3 and described in Table 3. Accordingly, the measured values per patient and the mean and standard deviation of each morphological measurement across all patients are listed in Table 4. Note that even though the standard deviation of each measurement may seem rather small, such small changes can lead to significant variations in the overall shape of the jaw, indicating the availability of high geometrical variations among different jaws under the study.
3.5. The pipeline for generating FE model of a human jaw The flexibility of the FE method allows it to use a wide range of spatial discretizations [34]. We opt for an unstructured tetrahedral mesh as can be robustly generated using automatic meshing tools [18] and can lead to similar accuracy and running time when using high order elements as structured meshes [35].
To eliminate the manual geometry processing, quality meshing, or mesh decimation steps, we propose to directly use the exported surface meshes as an input to our pipeline. Figure 4 shows a comparison between the conventional labor-intensive approach and our method for developing FE models of a human jaw. For the reproducibility purpose, the pipeline is implemented based on the free open-access geometry processing library libigl [19] and the meshing algorithm fTetWild [18]. The pipeline generates conformal volumetric meshes using imperfect meshes exported from the segmentation software with minimal human intervention. Figure 5 shows the flowchart of the utilized pipeline and the characteristics of the meshes at different steps.

Preprocessing
As can be seen in Figure 5, the input meshes to the pipeline are dense irregular meshes, which are not necessarily guaranteed to be watertight, manifold, and selfintersection-free, referred to as "triangle soup" in the computer graphics [23,22,36]. Hence, we apply a preprocessing step to both reduce the mesh sizes and produce meshes that are manifold, watertight, and self-intersection-free. These mesh characteristics assure meaningful values for the utilized signed distance functions in the next geometry processing steps of the used pipeline, i.e., the gap and PDL rim generation steps [36]. We use fTetWild [18] as a robust meshing tool to decimate and "clean up" the imperfect meshes. To be more specific, the teeth and bone geometries are tetrahedralized separately, and the boundary faces of the resulting tetrahedral meshes are then extracted as the cleaned-up reduced surface meshes to be used as the input meshes to the gap generation step. Further details on the surface mesh extraction process can be found in Section 3.5.6.

Gap generation for the PDL tissue
PDL has an average width of 0.2 mm [37]; its width can approximately be 0.15 mm around the middle third of the root and about 0.21 mm [38,39] to 0.38 mm [37] near the root apex and cervical regions.
Reconstructing the PDL layer using CBCT scans obtained in vivo from patients in clinics is a challenging process [40] as the commonly used voxel dimensions for the CBCT scans, ranging from 0.2 mm to 0.5 mm [3], are not fine enough to capture such a thin structure (roughly 0.2 mm) [41,42]. Although the geometry of the PDL layer can be reconstructed by segmenting it from the micro-CTs acquired in vitro, the x-ray exposure in such scans is extremely high and harmful for the human body, and it is usually obtained from dead specimens. Therefore, we first conduct the teeth-bone segmentation from the CBCT scans and then apply geometry processing techniques to create a gap between the teeth and bone geometries where the PDL can reside with an average thickness of 0.2 mm.
Ideally, to generate the required gaps for PDL geometries, we shrink the bone and teeth, each by 0.1 mm. First, we shrink the bone geometry by 0.1 mm, by moving its mesh points in the reverse direction of per-vertex normal with a magnitude of 0.1 mm, which we call it explicit shrinking approach. Before performing any explicit shrinking process, the radius curvature is locally computed for bone vertices to identify the sharp and thin structures. This is done by calculating the radius r of the mean curvature h at point t based on r(t) =1/h(t), as the reciprocal of the curvature at that point. The radius of the curvature provides useful information about the maximum magnitude that the surface vertices can be moved in the opposite direction of the normals of the surface before a singularity occurs. This means that the movement of the nodes with larger values would cause self-intersection issues and artifacts on the surface mesh. For this reason, evaluating the shrinking limits prior to the bone shrinking process is important. In the case where the shrinking limit is less than 0.1 mm, the pipeline shrinks the bone to its maximum shrinking limit, while more shrinking the teeth to compensate for the total desired gap.
To shrink the teeth, we use an implicit shrinking approach based on signed distance functions presented in Equation (1). The signed distance field contains the boundary of the geometry, i.e., zero iso-surface [43] and the information from different offset surfaces with positive and negative values. The positive offsets/iso-surfaces represent a dilated version of the geometry, while the negative values represent the eroded geometry. Consequently, we use an iso-contour of φ Ω (x) = −0.1 to shrink the tooth with a magnitude of 0.1 mm. Since the iso-contour is still an implicit representation of the shrunk tooth, we use a contouring method based on a marching cubes algorithm Table 4: Assessment of the geometrical variations of the mandibles and maxillae based on the computed measurements illustrated in Figure 3. Note that even small changes in the values can change the overall shape of the jaw considerably, indicating that high geometrical variations are available within our models.  [44] to convert it to an explicit representation including the vertices and connectivity matrix. As a result, wherever the radius curvature might be less than the desired offset value, the implicit shrinking approach can prevent any singularities at the root apexes. It can therefore be used as a more robust shrinking approach, especially for thin and sharp geometries with lower radius curvature values (here, the root apex of incisors). Note that the applied approach to the teeth creates a slightly wider space than the desired gap in the root apexes, which is in line with the clinical studies [45,39].

Boundary representation of PDL
Instead of an explicit representation of PDL which uses closed surface meshes, we use a boundary representation (B-Rep) approach using triangle meshes to describe the PDL domain. This needs to be distinguished from the basis spline (B-Spline) representation. we use B-Rep for describing the PDL domain for three main reasons. First, fTetwild uses the winding number information together with the input surface mesh to generate a volumetric mesh with no prior assumptions such as watertight or closed surfaces. Therefore, a shell surface can be used as a part of the input mesh to represent the domain boundary. Second, since fTetwild uses an implicit algorithm, it enables us to use B-Rep and perform boolean operations on different components to describe the domain boundary [18]. This is while it is not trivial to correctly define the B-Rep models for complex geometries using the Delaunay-based algorithms such as TetGen [46]. Third, B-Rep models help us to achieve zero-congruent contacting interfaces and avoid numerically small values for ε produced due to machine/floating-point precision. This, in turn, assures congruency and mesh conformity at the contacting interfaces for modeling the PDL layer as mentioned in Section 2.3.
The PDL geometry can be described using three main surfaces: the top surface of PDL, tooth-PDL, and PDLbone interfaces, in which the two last interfaces can be replaced with the teeth and bone geometries. To represent the free surface of the PDLs that are not in contact with the teeth and bone geometries, we generate only the top surface of each PDL, called the PDL rim in this study. B-Rep model using PDL rim and teeth-bone geometries as inputs to fTetWild helps to generate multi-domain volumetric mesh using boolean operations to guarantee geometry congruency and mesh conformity in the contacting interfaces.
We utilize a gap filling method [47] to generate the rims that connect the bone to shrunk teeth. This method detects an initial surface on the bone within the desired distance like 0.3 mm. Next, the boundary nodes of the detected surface are smoothened using a Laplacian smoothing approach to provide a base surface on the tooth socket with a smooth boundary on the bone. The smoothened boundary points are then snapped to the bone surface to assure the smoothened boundary is on the bone. Afterward, the base function with smoothened boundary is extruded with Figure 4: Finite element models created based on irregular dense meshes (left subfigure, A) exported from the segmentation step. Right subfigure: A comparison between the results of the conventional meshing approach (top-row) and utilized pipeline (bottom-row). The conventional approach involves time-consuming and labor-intensive geometry and mesh processing tasks (B); this results in non-congruent contacting interfaces (F), and non-conformal meshes (G). Note that the spotted marks in F indicate the undesired gaps/penetrations in the contacting interfaces. In contrast, the utilized method generates multi-domain volumetric meshes directly using the input irregular meshes from the generated PDL rims (shown in red in C) and guarantees the interface congruency as well as the mesh conformity as depicted in H.
the desired thickness (0.2 mm), and the extruded surface is snapped to the tooth surface to assure the points are located on the tooth surface. Finally, the PDL rim is constructed by connecting the boundary points of the base and extruded surfaces. We only use the PDL rims without the base/extruded surface to avoid any redundant representations of tooth-PDL and PDL-bone interfaces having potential numerical errors. These errors can cause issues in the volumetric mesh generation process of fTetwild, when the value of the user-defined maximum deviation parameter called envelope is less than or equal to the produced error in the duplicated surfaces.

Volumetric mesh generation
A unified volumetric mesh is generated using the surface meshes of the teeth, bone, and produced PDL rims, It is obtained by applying a union operation on the provided input surface meshes using fTetWild [18] with optimal envelope ( ) and ideal edge length values of 2 × 10 −4 and 0.01, respectively. Applying a union operation using fTetWild generates a single volumetric mesh for the bounding box surrounding all input surface meshes called the raw mesh [18]. Note that the default filtering method in fTetWild exports a labeled version of the raw mesh that assigns the tetrahedra outside the closed input surfaces as the background elements.

Proposed tetrahedra filtering method
The default tetrahedra filtering/labeling method in fTetWild can only be used when all input meshes are closed surface meshes, which is not the case for the PDL rims. Hence, a specific filtering approach is used for as-signing each tetrahedron a label associated with one of the input domains, i.e., the teeth, PDL, or bone. The remained tetrahedra are labeled as background. We use a distance field-based labeling method for teeth and bone. To do so, the barycenter of each tetrahedron is used as an average point of tetrahedron vertices to decide whether the element is positioned inside or outside a surface considering the distance field of each barycenter with respect to each of the teeth and bone geometries. The negative distance values indicate the points are located inside the surface mesh.
We propose using a signed distance-based method combined with the winding number information to obtain the PDL geometries. The PDL is expected to reside in the gap between the teeth and bone; hence, as shown in Figure 6, we first select all tetrahedra with positive distance values below 0.3 mm around teeth (A) and bone (B) based on the intersection of these two sets of tetrahedra (C). Next, the intersection of the obtained tetrahedra (C) and the tetrahedra with positive winding numbers with respect to the PDL rims (D) is utilized to produce a smooth boundary (E) for the top-free surface of each PDL. The winding numbers, as also used in fTetWild for implicit meshing, are applied to eliminate all meshing artifacts introduced by the intersection operation (C). As can be seen, applying the proposed filtering method considering the information of the winding numbers to the raw mesh provides smooth PDL top surfaces as well as conformal meshes in the tooth-PDL and PDL-bone interfaces, which include shared nodes in the contacting interfaces and provide perfect adhesions in the contacting interfaces. the contacting interfaces and provide perfect adhesions in the contacting interfaces.

Surface mesh extraction
After labeling all tetrahedra in the raw mesh, we select boundary faces of each domain, i.e., the shared faces between the domain and background tetrahedra, and save them as reduced adaptive surface meshes. Note that the extracted surfaces are zero-congruent and conformal meshes with no meshing artifacts.

FE simulation setup
The used pipeline automatically sets up the FE problems based on two different scenarios; (1) uncontrolled tipping and (2) biting scenarios. It automatically generates FEM files of the uncontrolled tipping scenario for all patients and FEM files of the biting scenario for those whose scans were acquired in a natural biting position. Note that the utilized code for defining the boundary and loading conditions is scenario-specific, which requires to be adapted for other scenarios. The defined FE problem for each scenario is exported as an XML file including the tetrahedral meshes, types of elements, utilized material models and properties of each domain, as well as boundary and loading conditions.
In these scenarios, we use Tet4 elements for the teeth and bone with negligible deformations. Furthermore, to avoid the locking artifacts [48] that are common when linear tetrahedral elements are used to model large deformations, we apply a quadratic basis to the elements of the PDL layer, to consider Tet10 elements for the PDL domain. We prefer to use an explicit tetrahedral mesh to model the PDL layer [49,40,16] (instead of the cheaper option of using shell elements) as it can more faithfully capture shear, bending, and buckling of the membrane, and at the same time, it simplifies both meshing and simulation, as it does not require special handling for inserting shell elements in a tetrahedral mesh and coupling between the volumetric elastic model and its interaction with the shell. Further, providing a volumetric representation of the PDL layer allows future studies to explore the effects of utilizing visco-elastic-plastic material models for the PDL layer.

Surface mesh extraction
After labeling all tetrahedra in the raw mesh, we select boundary faces of each domain, i.e., the shared faces between the domain and background tetrahedra, and save them as reduced adaptive surface meshes. Note that the extracted surfaces are zero-congruent and conformal meshes with no meshing artifacts.

FE simulation setup
The used pipeline automatically sets up the FE problems based on two different scenarios; (1) uncontrolled tipping and (2) biting scenarios. It automatically generates FEM files of the uncontrolled tipping scenario for all patients and FEM files of the biting scenario for those whose scans were acquired in a natural biting position. Note that the utilized code for defining the boundary and loading conditions is scenario-specific, which requires to be adapted for other scenarios. The defined FE problem for each scenario is exported as an XML file including the tetrahedral meshes, types of elements, utilized material models and properties of each domain, as well as boundary and loading conditions.
In these scenarios, we use Tet4 elements for the teeth and bone with negligible deformations. Furthermore, to avoid the locking artifacts [48] that are common when linear tetrahedral elements are used to model large deformations, we apply a quadratic basis to the elements of the PDL layer, to consider Tet10 elements for the PDL domain. We prefer to use an explicit tetrahedral mesh to model the PDL layer [49,40,16] (instead of the cheaper option of using shell elements) as it can more faithfully capture shear, bending, and buckling of the membrane, and at the same time, it simplifies both meshing and simulation, as it does not require special handling for inserting shell elements in a tetrahedral mesh and coupling between the volumetric elastic model and its interaction with the shell. Further, providing a volumetric representation of the PDL layer allows future studies to explore the effects of utilizing visco-elastic-plastic material models for the PDL layer. Material models and parameters. In this study, the bone and teeth tissues deform negligibly under the applied forces. Therefore, we assume no distinctions between different structures of the bone (cortical and trabecular) and teeth (enamel, dentin, and pulp) [4,50,51]. Besides, the porous fibrous periodontal ligament tissue is assumed as a homogeneous structure [49]. We use the Neo-Hookean material model with Poisson's ratios of 0.3, 0.3, and 0.45, and Young's modulus values of 2000 MPa, 1500 MPa, and 68.9 MPa to describe the mechanical behavior of the tooth, PDL, and bone domains, respectively [4,52,2].
Note that the used simplex material models can be replaced with any complex constitutive models to properly mimic the anisotropic viscoelastic behavior of the PDL [51,53,49], or orthotropic characteristic of the bone. This can however increase the computational costs of the simulations.
Boundary conditions. In the tipping scenario, a Dirichlet boundary condition is defined on all nodes located at the bottom/top surface of the mandible/maxilla to fix displacements of the nodes in all three directions. In the biting scenario, the same Dirichlet boundary condition is applied on the top nodes of the maxilla to fix them in all three directions. Besides, a similar Dirichlet boundary condition is used to restrict the movement of the mandible only in the anterior-posterior and medial-lateral directions with imposing no restrictions in the third direction.
Loading conditions. In the tipping scenario, we apply a perpendicular force with a magnitude of 1 N at the center of each crown to mimic the uncontrolled tipping motion in the lingual direction. In the biting scenario, a pressure load of 200 N is applied to the bottom surface of the mandible to simulate the biting force.

Contact definition.
To have a perfect adhesion in the tooth-PDL and PDL-bone interfaces, we generate a single volumetric mesh by combining surface meshes of the shrunk teeth, PDL rim, and bone geometries using a union operation in fTetWild [18]. The nodes at the contacting interfaces are shared between the adjacent domains, thus guaranteeing a complete adhesion as well as conformal meshes at the contacting interfaces. Therefore, there will be no sliding, separation, or penetration at the tooth-PDL and PDL-bone contacting interfaces. Further details on the volumetric mesh generation can be found in Section 3.5.4. Furthermore, in the both scenarios, any potential contacts between different teeth are modeled using the incremental potential contact formulation [54].

Technical details
Implicit shrinking approach. The explicit representation of geometries only provides boundary information of them. This is while the implicit representation of geometries based on the signed distance function provides useful information about the interior and exterior of the shape as well as the outline/boundary of geometry. Therefore, we use boundary information obtained from the explicit representation to define the signed distance function/filed before evaluating it on the desired points in an R 3 space and shrinking the teeth. To do so, we sample the distance field in the R 3 space using a regular 3D grid in the bounding box of each tooth. A fine grid sampling size of 0.1 mm is used to be smaller than the minimum isotropic voxel size (0.15 mm) as presented in Table 2, to avoid the aliasing effect or spiky surfaces in the sampling process and have smooth geometries in the shrunk teeth.
Next, the signed distance values are obtained for each point of the grid considering the surface mesh of the tooth. Note that sampling using regular grids requires excessive memory for geometries with large dimensions, which is not the case here as the dimension of each tooth is relatively small. In general, for meshes with fine details that cover a larger space, it is suggested to use adaptive sampling approaches such as the octree-based methods [55], to densely sample voxels in specific regions near the boundary or regions with great details.
Alternative PDL modeling. An alternative approach for generating FE models of the tooth-supporting complex can be using Delaunay-based volumetric meshing tools like TetGen. Still, one needs to assure mesh conformity or interface congruency on the surface meshes, and next, to enforce the explicit volumetric meshing algorithm to preserve the surface meshes [10]. Note that providing quality surface meshes with mesh conformity criteria per se is not trivial. Besides, even though the alternative explicit approach can provide quality surface meshes, it may not provide optimal quality for the tetrahedra due to the additional constraint applied to preserve the surface mesh on the boundary of the domain.
Advantages of implicit over explicit meshing algorithms. The explicit volumetric mesh generation approaches such as the Delaunay-based algorithms with incremental local mesh operations generate tetrahedral meshes covering the domain interior and are highly faithful to the input mesh. This means that the algorithms cannot coarsen the triangles on the input surface meshes. Therefore, they cannot produce coarse quality tetrahedra on the boundary faces of the input mesh, indicating that the element size at the domain boundary depends on the dense and irregular/rugged surfaces exported from the segmentation step. Consequently, one needs to prepare an adaptive surface mesh with high-quality triangles before applying the Delaunaybased algorithms along with a sizing field function to cre- ate adaptive volumetric mesh both in the domain interior and its boundary. Furthermore, since the explicit meshing algorithms fail in generating volumetric meshes in the presence of any self-intersections, one needs to clean up the meshes before applying the Delaunay-based meshing algorithms to create quality and adaptive mesh. In contrast, implicit volumetric meshing approaches like fTetWild impose no assumptions on the input mesh and can handle imperfect input meshes with selfintersection artifacts. They can also produce an adaptive mesh that provides coarse tetrahedra at the domain boundary, slightly deviating from the input surface mesh, based on a user-defined input parameter. In addition, fTetWild uses winding numbers along with surface meshes to generate tetrahedral meshes based on an implicit meshing approach. This enables us to use the constructive solid geometry (CSG) model for describing a domain with complex boundaries by combining several simpler domains using boolean operations [56]. This is while the explicit meshing tools cannot support CSG models [46].
Volumetric meshing. The fTetWild algorithm uses the user-defined epsilon and ideal edge length parameters to control the output mesh's accuracy and size. The epsilon value indicates how much fTetWild can deviate from the input surface mesh for generating a quality volumetric mesh. Hence, using a smaller epsilon value preserves the input geometry's details at the cost of higher computational time. In contrast, larger values can provide less accurate meshes in less computational time. Moreover, smaller values of ideal edge length provide denser meshes, Figure 7: Comparison of the quality histograms of four different quality measurements used for the evaluation of the quality of the tetrahedra generated by OpenMandible and Open-Full-Jaw. Left to right: The radius-edge ratio Q rl ; The volume-edge ratio Q vl ; The radius ratio Q rr ; The angle measurement Q θ . Note that L max , R in , and R out denote the maximum edge length, the radius of the insphere, and radius of the circumsphere of a tetrahedron, respectively, and V and S min represent the volume and minimum face area of the tetrahedron.
these two parameters can control the output mesh size and its geometrical accuracy. In this work, the utilized values for epsilon and ideal edge length are obtained based on a grid-search method in the intervals [10 −6 , 10 −3 ] and [0.005, 0.02], respectively. Note that the fTetWild multiplies these parameters by the length of the diagonal of the input meshes' bounding box. Hence, the parameters are sensitive to the size of the bounding box that includes all input meshes. In other words, they need to be adjusted according to the model's size; if, for instance, the model is a cut-out model with a much smaller bounding box dimension than the full jaw model.

Results and discussion
In this section, the generated models are assessed from different aspects, and the results are compared with the state-of-the-art. In Section 4.1, we measure the thickness of the generated PDL layers. In Section 4.2, different mesh properties such as mesh sizes and mesh qualities are evaluated for the developed models; the results are compared to the OpenMandible model. In Section 4.3, FE simulation results of different patients are presented under identical tipping and biting scenarios.

Generated PDL properties
We assume an average thickness of 0.2 mm for the PDL layers according to the literature [7,9] and assess the thicknesses of the generated PDLs across different patients models by computing the distances between the points located on the outer and inner surfaces of the PDL geometries. The thickness details of the generated PDLs can be seen in Table 5. As can be noticed from the table, the computed minimum, maximum, and average PDL thickness values are in line with those reported in the literature, i.e., 0.15, 0.3, and 0.2mm, respectively [37].
Likewise, we investigate the thicknesses for the PDLs from the OpenMandible dataset. As it can be seen in the last row of Table 5, the obtained thicknesses for Open-Mandible are about two-three times the values reported in the literature. As the PDL thickness plays an essential role in exerting the applied load from teeth surfaces to the adjacent bone and teeth movements, it can be deduced that under identical scenarios and loading systems the Open-Mandible model cannot result in stress/strain values comparable with the ones obtained in this study.

Mesh properties
We evaluate the total number of elements and mesh qualities of the developed models and compare the mesh qualities of our model with those of the OpenMandible dataset. The maximum mesh deviation between the input and output meshes of ftetwild is then calculated to examine whether the meshing process alters the geometries negligibly.

Mesh sizes
We develop all computational meshes of our study on a machine with a 2.60 GHz processor and 16 GB of RAM, which takes around 40.78 ± 12.19 minutes per jaw. The details of the original dense meshes, the generated volumetric meshes, and their extracted surface meshes are summarized in Table 6. The pipeline produces conformal adaptive meshes from irregular dense meshes. The utilized adaptive meshing approach helps to reduce the total number of the elements, by producing coarser elements in regions far from the teeth while generating finer elements on root apexes, PDL layer, alveolar crest, or regions with fine details like thin walls of the maxilla.

Mesh quality analysis
The mesh quality required for an FE analysis can vary depending on the application and utilized numerical meth-14 Figure 7: Comparison of the quality histograms of four different quality measurements used for the evaluation of the quality of the tetrahedra generated by OpenMandible and Open-Full-Jaw. Left to right: The radius-edge ratio Q rl ; The volume-edge ratio Q vl ; The radius ratio Q rr ; The angle measurement Q θ . Note that L max , R in , and R out denote the maximum edge length, the radius of the insphere, and radius of the circumsphere of a tetrahedron, respectively, and V and S min represent the volume and minimum face area of the tetrahedron.
while larger values can produce coarser elements. Hence, these two parameters can control the output mesh size and its geometrical accuracy. In this work, the utilized values for epsilon and ideal edge length are obtained based on a grid-search method in the intervals [10 −6 , 10 −3 ] and [0.005, 0.02], respectively. Note that the fTetWild multiplies these parameters by the length of the diagonal of the input meshes' bounding box. Hence, the parameters are sensitive to the size of the bounding box that includes all input meshes. In other words, they need to be adjusted according to the model's size; if, for instance, the model is a cut-out model with a much smaller bounding box dimension than the full jaw model.

Results and discussion
In this section, the generated models are assessed from different aspects, and the results are compared with the state-of-the-art. In Section 4.1, we measure the thickness of the generated PDL layers. In Section 4.2, different mesh properties such as mesh sizes and mesh qualities are evaluated for the developed models; the results are compared to the OpenMandible model. In Section 4.3, FE simulation results of different patients are presented under identical tipping and biting scenarios.

Generated PDL properties
We assume an average thickness of 0.2 mm for the PDL layers according to the literature [7,9] and assess the thicknesses of the generated PDLs across different patients models by computing the distances between the points located on the outer and inner surfaces of the PDL geometries. The thickness details of the generated PDLs can be seen in Table 5. As can be noticed from the table, the computed minimum, maximum, and average PDL thickness values are in line with those reported in the literature, i.e., 0.15, 0.3, and 0.2mm, respectively [37].
Likewise, we investigate the thicknesses for the PDLs from the OpenMandible dataset. As it can be seen in the last row of Table 5, the obtained thicknesses for Open-Mandible are about two-three times the values reported in the literature. As the PDL thickness plays an essential role in exerting the applied load from teeth surfaces to the adjacent bone and teeth movements, it can be deduced that under identical scenarios and loading systems the Open-Mandible model cannot result in stress/strain values comparable with the ones obtained in this study.

Mesh properties
We evaluate the total number of elements and mesh qualities of the developed models and compare the mesh qualities of our model with those of the OpenMandible dataset. The maximum mesh deviation between the input and output meshes of ftetwild is then calculated to examine whether the meshing process alters the geometries negligibly.

Mesh sizes
We develop all computational meshes of our study on a machine with a 2.60 GHz processor and 16 GB of RAM, which takes around 40.78 ± 12.19 minutes per jaw. The details of the original dense meshes, the generated volumetric meshes, and their extracted surface meshes are summarized in Table 6. The pipeline produces conformal adaptive meshes from irregular dense meshes. The utilized adaptive meshing approach helps to reduce the total number of the elements, by producing coarser elements in regions far from the teeth while generating finer elements on root apexes, PDL layer, alveolar crest, or regions with fine details like thin walls of the maxilla.

Mesh quality analysis
The mesh quality required for an FE analysis can vary depending on the application and utilized numerical meth- Table 6: The total number of elements (×10 3 ) of the surface meshes (triangles) and volumetric meshes (tetrahedra). The number of elements is significantly reduced in the adaptive output surface meshes compared to that of the irregular dense input surface meshes.
ods [57]. In general, a regular tetrahedron has the highest mesh quality in computational models, and the main guideline is to avoid using low-quality/badly-shaped tetrahedra, as they can affect the accuracy of the numerical methods.
The OpenMandible uses TetGen to obtain volumetric meshes from the manually generated conformal surface meshes. It sets the upper limit of the radius-edge ratio of to-be-generated tetrahedra to 1.5. This mesh quality constraint controls the ratio between the radius of the circumscribed sphere and the shortest edge of each tetrahedron, to prevent the production of low-quality/badly-shaped tetrahedra. In addition to this quality constraint, Open-Mandible enforces TetGen to preserve the provided input surface meshes to have conformal volumetric meshes. Preserving the surface mesh is the only approach to produce conformal volumetric meshes when using explicit volumetric meshing approaches. This raises the question of whether TetGen can achieve the specified quality constraint value while enforcing another restriction to preserve the input surface mesh.
We quantitatively assess the quality of the generated volumetric meshes of this study and those of the Open-Mandible. To be more specific, four different quality measurements presented in [57], i.e., the radius-edge ratio Q rl [58], the volume-edge ratio Q vl [59,60], the radius ratio Q rr [61,62], and the angle measurement Q θ [61], are used to evaluate the quality of the volumetric meshes. For a fair comparison, we compare the results from one of our patients to those of the OpenMandible study, as shown in Figure 7.
In a regular tetrahedron, the values of each utilized quality measurement correspond to one, indicating that an ideal high-quality mesh is expected to have a histogram peak at one. Therefore, the resulting distributions (normal or skewed normal) of all quality histograms show that our generated volumetric mesh has higher quality elements compared to OpenMandible. The proposed model also has narrower distributions with small discrepancies around their means. This holds even in the last case (Q θ ), where one of the peaks of the distribution (mode) for OpenMandible is closer to one. Moreover, the Open-Mandible model sees two or more peaks in its distributions that can be modeled by mixed normal distributions.

Mesh deviation analysis
If one does not enforce the surface preservation criterion, the Delaunay-based volumetric meshing algorithms like TetGen include all points of the input surface mesh and a number of additional points. Hence, the generated volumetric meshes have boundary elements as fine as the input surface mesh. In contrast, fTetwild provides volumetric meshes as coarse as possible on the surface while preserving the input geometry. This algorithm slightly deviates from the input surface according to the user-defined maximal deviation (envelope) value. Therefore, to generate an accurate volumetric mesh, we must not deviate too much from the input surface mesh.
Accordingly, we quantitatively evaluate the deviation of final extracted surface meshes from the input surface meshes. We use Hausdorff distance as an error measurement between the input and extracted surface meshes from the generated volumetric meshes. The obtained distance values are presented in Table 7. As it can be seen, the maximum deviation of the meshes both in mandibular and maxillary jaws and their associated teeth is negligible with respect to the dimension of the entire jaw models.

Mesh deviation analysis
If one does not enforce the surface preservation criterion, the Delaunay-based volumetric meshing algorithms like TetGen include all points of the input surface mesh and a number of additional points. Hence, the generated volumetric meshes have boundary elements as fine as the input surface mesh. In contrast, fTetwild provides volumetric meshes as coarse as possible on the surface while preserving the input geometry. This algorithm slightly deviates from the input surface according to the user-defined maximal deviation (envelope) value. Therefore, to generate an accurate volumetric mesh, we must not deviate too much from the input surface mesh.
Accordingly, we quantitatively evaluate the deviation of final extracted surface meshes from the input surface meshes. We use Hausdorff distance as an error measurement between the input and extracted surface meshes from the generated volumetric meshes. The obtained distance values are presented in Table 7. As it can be seen, the maximum deviation of the meshes both in mandibular and maxillary jaws and their associated teeth is negligible with respect to the dimension of the entire jaw models.
16 Figure 9: The obtained displacement (A) and stress (B) fields for the mandibular jaw of Patient 17 under the tipping scenario. Note that stress concentration is observed on the lingual side of anterior teeth due to the labiolingual tipping load.

Mesh deviation analysis
If one does not enforce the surface preservation criterion, the Delaunay-based volumetric meshing algorithms like TetGen include all points of the input surface mesh and a number of additional points. Hence, the generated volumetric meshes have boundary elements as fine as the input surface mesh. In contrast, fTetwild provides volumetric meshes as coarse as possible on the surface while preserving the input geometry. This algorithm slightly deviates from the input surface according to the user-defined maximal deviation (envelope) value. Therefore, to gener-ate an accurate volumetric mesh, we must not deviate too much from the input surface mesh.
Accordingly, we quantitatively evaluate the deviation of final extracted surface meshes from the input surface meshes. We use Hausdorff distance as an error measurement between the input and extracted surface meshes from the generated volumetric meshes. The obtained distance values are presented in Table 7. As it can be seen, the maximum deviation of the meshes both in mandibular and maxillary jaws and their associated teeth is negligible with respect to the dimension of the entire jaw models.

FEM verification
As the last part of our analyses, we assess the displacement and stress fields in the tipping and biting scenarios to ensure that the stress patterns are smooth and have no unrealistic stress concentrations. To do so, we run the simulations using the PolyFEM [63] FE solver. PolyFEM is an FE simulation toolkit that supports elastodynamic deformations with linear and non-linear material models. It provides an adaptive p-refinement that allows increasing the order of basis functions for specific domains while utilizing linear basis functions for the other domains. Hence, we use Tet10 elements for the PDL layer to increase the simulations' accuracy and avoid element locking issues [64]. Besides, PolyFEM uses the incremental potential contact formulation [54] for contact response and friction, which ensures valid, penetration-free meshes during the entire simulation. The contacts are automatically detected by proximity; hence, there is no need to specify contact surfaces, which significantly simplifies the scene setup.
The simulation results are visualized in ParaView [65]. Figure 8 shows the result of the biting scenario of a selected patient, including the resulting displacement and stress fields. Due to the deformation of the PDL layer under the applied biting load, the displacement fields can be seen both on the mandibular jaw and the posterior maxillary teeth (A). Besides, stress concentrations can be observed on the associated PDL layers, especially at root apexes and bifurcation regions indicated by arrows in (B). Figure 9 illustrates the result of the auto-generated uncontrolled tipping scenario for the mandibular jaw of the same patient. As seen in [4], the anterior teeth have higher displacement fields than the posterior teeth under an identical load (A). Hence, higher stress concentrations can be seen on the lingual side of the PDLs associated with the anterior teeth (B). Additionally, simulation results of different patients in the tipping scenario are shown in Figure 10. As can be noticed, the stress values considerably change from one patient to another, indicating the importance of utilizing population models for multi-patient analysis.
It should be noted that although we have tested our models under two scenarios, the developed volumetric meshes can be used in various scenarios and different FE frameworks. Furthermore, the FE models can benefit from using more complex material models, boundary, and loading conditions. For example, the provided meshes can be integrated with outputs of other studies [10,9] to consider masticatory muscles for more realistic biting scenarios.

General discussion
The utilized and conventional meshing approaches generate the volumetric meshes using reconstructed geometries based on accurately segmented scans. However, obtaining such an accurate segmentation is inherently timeconsuming and labor-intensive and, in some cases, could be highly challenging due to the complexity of the problem or lack of high-resolution scans [66]. The templatebased deformation techniques [67,68] can be used to automatically reconstruct the 3D geometries by creating a template mesh and deforming it according to the new data samples. Still, deforming a template mesh using registration approaches or deep learning methods requires accurate spatial registration and high-quality volumetric meshes with no distorted elements.
Obtaining accurate spatial registration and high-quality volumetric meshes for the human jaw can be challenging, as there are large variations among geometries of different patients (Figure 2), such as geometrical differences in the bones and teeth, missing teeth, and topological changes in the number of roots, e.g., the mandibular and maxillary molars with two to four roots. Therefore, to include different types of variations in the data for a plausible deformation, one needs to have different templates covering missing teeth or various numbers of roots. This in turn is a time-consuming process and can increase the complexity of the model. On the other hand, large deformations in small volumes of teeth and roots can result in distorted elements, preventing us from generating highquality meshes, especially in the PDL layers with thin structures that need to be modeled with fine volumetric elements.
The current study introduces the largest-ever dataset of patient-specific human jaws reconstructed from CBCT scans. We believe this unique clinically validated dataset would pave the way for future population studies in the field. More specifically, data augmentation techniques using machine learning [69,70] can be applied to the Open-Full-Jaw dataset to expand its size and variability by generating plausible synthetic data. In addition, this would enable us to use deep learning methods, which require a large amount of data for training. Still, one needs to use a dataset with enough variations for sampling and assess the generated samples' validity.   Figure 10: A comparison of stress distributions of all utilized jaws in the tipping scenario with identical load magnitudes. The stress fields are clipped within the range [0, 1] MPa to represent the changes better using the same color map. As can be noticed, the stress values and concentrations considerably vary across different patients, which indicates the importance of utilizing population models for multi-patient analysis and generalizability in computational modeling studies.

Conclusion
18 Figure 10: A comparison of stress distributions of all utilized jaws in the tipping scenario with identical load magnitudes. The stress fields are clipped within the range [0, 1] MPa to represent the changes better using the same color map. As can be noticed, the stress values and concentrations considerably vary across different patients, which indicates the importance of utilizing population models for multi-patient analysis and generalizability in computational modeling studies.
diku-dk/Open-Full-Jaw), with patient-specific models of 17 human mandibles and maxillae. The dataset contains clinically validated segmented geometries shared as dense surface meshes and adaptive quality volumetric with conformal meshes in the contacting interfaces. It also includes the principal axes for each patient's teeth and the generated FEM files of the uncontrolled tipping and biting scenarios for all patients. Finally, we share the nearly-automated pipeline used for geometry processing, re-meshing, and generating volumetric meshes.
In addition, we evaluated the generated models and quantified them in terms of the mesh quality and accuracy of the models, and compared the results with the stateof-the-art. The obtained results indicate that the developed computational models are precise, considering the low error/distance from input surface meshes. Moreover, the quality of the volumetric elements evaluated based on different quality measurements imply that the generated volumetric meshes consist of quality elements suitable for the FEM of the human jaw. Hence, we believe the Open-Full-Jaw dataset can be used in various FE scenarios and a wide range of intra-and inter-patient analyses.
The shared repository includes all detailed information for reproducing the models of this study. In addition, the utilized pipeline allows other researchers in the field to generate quality volumetric meshes and FE model files directly using dense and irregular meshes with minimal human intervention. This will help other researchers easily extend their datasets without spending much time and effort on manually cleaning up the meshes and non-trivially producing conformal meshes. Furthermore, similar concepts as those used in this study to generate population models of the tooth-supporting complex can be adapted to other areas, such as pelvic girdles and hip joints [71].

Funding
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 764644. This paper only contains the authors' views, and the Research Executive Agency and the Commission are not responsible for any use that may be made of the information it contains. 3Shape A/S provided financial support in the form of salaries for authors TG and PS. The funders had no role in study design, data collection/analysis, publication decision, or manuscript preparation. This work was also partially supported by the NSF CA-REER award under Grant No. 1652515, the NSF grants OAC-1835712, OIA-1937043, CHS-1908767, CHS-1901091, NSERC DGECR-2021-00461 and RGPIN-2021-03707, a Sloan Fellowship, a gift from Adobe Research and a gift from Advanced Micro Devices, Inc.

Ethical approval
All data used in this paper was provided, already anonymized, by 3Shape A/S. The data was originally acquired for diagnostic purposes unrelated to this study. No other aspect of this work triggered ethical issues.