2. Materials and Methods
The study was conducted through a collaboration between the Maxillofacial Surgery Unit of AOU Città della Salute e della Scienza, University of Turin, Italy, and the Department of Management and Production Engineering (DIGEP), Polytechnic University of Turin, from October 2025 to March 2026.
Patients’ selection
Dental models from patients with occlusal discrepancies and undergoing orthognathic surgical treatment between October 2025 and January 2026 were collected. For each patient, preoperative dental arches and surgical occlusion defined using the manual method on physical plaster casts (Manual Reference Occlusion, MrO), was acquired in digital format through intraoral scanning (IOS); the resulting three-dimensional (3D) meshes were exported in .STL format for further processing.
The inclusion criteria for the analysis were (1) age ≥ 18 years; (2) presence of an occlusal discrepancy requiring surgical correction; (3) availability of preoperative dental models in .STL format. Exclusion criteria were applied as follows: (1) noisy or incomplete IOS acquisitions; (2) presence of fewer than six pairs of antagonistic teeth per side; (3) segmental surgery of the maxilla; (4) inability to retrieve clinical information.
For each subject, demographic (gender and age) and clinical data (type of facial dysmorphism and surgical procedure) were collected, together with the .STL file representing the dental arches and final MrO, which was considered the reference standard. The dataset comprised patients with a broad spectrum of skeletal and dental malocclusions. No subgroup analysis based on the surgical approach (single jaw versus bimaxillary surgery) was conducted, as the definition of the final occlusion is independent of the number of jaws undergoing surgical repositioning.
Based on the data collected from these patients, a digital surgeon-oriented interactive interface was developed incorporating a two-step workflow: 1) a geometry-based pipeline for initial automatic positioning of the dental arches followed by 2) a dedicated Graphic User Interface (GUI) enabling the operator to adjust the occlusion through a surgeon-oriented standardized refinement phase.
Geometry-based automatic pipeline
The first step of the proposed methodology relies on the development of an automatic pipeline for the estimation of the best occlusal relationships between dental arches based on geometric references on 3D models.
Firstly, 3D meshes from IOS were processed and optimized using MeshLab (version 2023.12; ISTI-CNR, Italy). Subsequent computational algorithms for alignment and quantitative analysis were developed in MATLAB (version R2023a; MathWorks, Natick, MA, USA), allowing full control over computational procedures and ensuring reproducible results.
The high level of geometric complexity introduced by detailed IOS makes direct application of computational algorithms challenging. To reduce geometrical complexity, the dental surfaces were represented in form of depth maps preserving the morphological information necessary to assess inter-arch contacts and relationships without reproducing the entire dental geometry, reducing data dimensionality and sensitivity to non-functional variations.
To compute the depth maps, each mesh was sampled on a regular grid of 512x512 points and for each point the depth value was calculated using a volumetric ray casting algorithm [
32]. At this point, the 3D geometry was projected into a 2D domain while retaining relevant occlusal information, making subsequent processing more efficient [
33].
To isolate relevant geometry from areas of non-interest, the depth maps were cropped (
Figure 1).
These simplified models were standardized using principal component analysis (PCA), ensuring consistent orientation and alignment of the datasets in a common reference system [
34].
Curvature-based descriptors were then computed on these surfaces, specifically the mean curvature (H) and the Gaussian curvature (K), from which the principal curvatures are derived. These parameters ultimately enabled the calculation of the shape index (S), providing a local classification of surface morphology and allowing the identification of concave, convex, and saddle-shaped regions [
35].
To ensure consistency between the upper and lower arches, the sign of the mandibular shape index was inverted to account for the opposite spatial orientation of the mandibular arch (
Figure 2). At this point, the shape index values allow the identification of matching points between the two arches, by pairing the concavities of one arch with the convexities of the other, and vice versa.
The resulting maps were stored and used as input features for the subsequent alignment process.
The mandibular arch was considered as a fixed reference, while the maxillary arch was considered as a mobile model. An initial rigid transformation of the arches was obtained using a 2D variant of the Iterative Closest Point (ICP) algorithm, which uses a least-squares approach based on the Kabsch algorithm [
36,
37]. The algorithm analyses the points of the depth maps, identifies the corresponding points based on the shape index, and subsequently minimizes the distance between them, thereby progressively bringing the two arches closer together by translating them along the x- and y-axis while neglecting translation along the z-axis.
The adoption of the 2D ICP algorithm, optimally suited for shape information [
38], allows the reduction of computational complexity without losing the complex geometrical information of the anatomical models, to which 3D ICP is excessively sensitive. The 2D ICP applied to shape index allows to couple the cusps of the maxillary arch with the corresponding fossae of the mandibular arch providing a global reciprocal positioning between maxillary and mandibular models (
Figure 3).
To combine the functionality of geometric ICP with the specific functionality of the shape index, a weighted distance between corresponding points on the two arches is defined by combining the Euclidean distance of the ICP with that related to the shape index. This definition is formalized in the following formula:
where:
is the Euclidean distance on the XY plane;
is a similarity distance based on the shape index;
is the weight accounting for morphological consistency.
In the present study, an experimental α value of 0.3 was adopted, meaning that the shape index contributed 30% to the overall distance calculation; preliminary tests on a subset of cases demonstrated that α = 0.3 provided the best compromise between geometric proximity and morphological correspondence. By adjusting this parameter, the algorithm can range from a purely geometric ICP-based behaviour (α = 0), to a fully shape index–based behaviour (α = 1).
After applying the estimated transformations, a specific refinement phase was implemented to correct the translation on the Z-axis, corresponding to the vertical translation (Z-shift) required to achieve a physiologically consistent inter-arch distance.
To achieve this, a preliminary correction of the rotation around the X-axis was required in order to obtain a maxillary occlusal plane parallel to the mandibular one, thereby ensuring a uniform Z-axis distance between the corresponding points of the two arches. To do so, the algorithm divides the model along the anteroposterior axis into anterior (incisor) and posterior (molar) regions. Representative surface points on the maxillary arch are selected within these regions, and their vertical distances from the mandibular arch are evaluated. The height difference between the anterior and posterior regions is then used to estimate the corrective inclination, which is subsequently applied as a rotation around the X-axis to compensate for occlusal plane misalignment (
Figure 4).
Following tilt correction, the vertical translation (Z-shift) required to achieve a physiologically consistent inter-arch distance is calculated by identifying corresponding points between the two arches via spatial proximity and analyzing their vertical differences. Robust statistical measures such as percentiles are used to reduce the influence of outliers and irregularities on the local surface. The Z-axis translation is applied as the final step before the refinement procedure.
To address cases where excessive inter-arch penetration occurs, a non-penetration constraint was introduced. Considering that teeth can undergo slight deformations under load, an acceptable penetration range of 300µm to 500µm was established during the experimental phase [
39,
40,
41]. In this study, an penetration threshold of 0.3 mm was adopted, allowing a small, controlled degree of interpenetration between the arches. However, the framework has been designed to enable the practitioner to modify this value.
Two methods for managing the non-penetration constraint are proposed below:
The algorithm evaluates the spatial relationship between the maxillary and mandibular arches by measuring the inter-arch distances. Areas of overlap are identified as penetrations. Based on this analysis, the extent and distribution of penetrations are calculated to determine the vertical adjustment required to obtain the desired clearance between the arches. The rigid constraint is then applied through a rigid vertical translation of the maxilla model, ensuring a controlled and consistent correction of the inter-arch relationship. The final alignment is subsequently verified to confirm that the required clearance has been achieved. This constraint, as shown in
Figure 5, is effective in spacing the arches but provides few stable occlusal contacts.
- 2.
This method, introduced to increase the number of occlusal contacts, formulates the alignment as an iterative optimization problem, in which a rigid transformation of the maxilla is computed to maximize the quality of the occlusion.
The optimization is based on an objective function defined as a weighted combination of three terms:
a penalty for penetration below a target free space value;
a contact reward modeled by a Gaussian function centered on the free space value;
a regularization term to limit non-physiological transformations.
Moreover, a spatial weighting strategy was introduced to modulate the contribution of the three terms along the dental arch, assigning greater importance to the posterior regions, particularly the molars and premolars, while preserving a smaller, yet non-zero, contribution from the anterior region. This weighting is implemented through a logistic function governed by parameters defining the transition point, slope, and relative weighting between the anterior and posterior regions.
At each iteration, the maxillary model is rigidly transformed, and the point vertical spaces are recalculated with respect to the mandibular model. Optimization favors configurations that maintain stable contact while avoiding excessive separation or penetration, resulting in smooth “sliding” behavior between the arches.
The iterative process is performed for up to 900 iterations or until convergence is achieved, typically before the completion of the total number of iterations. This phase follows the initial alignment and allows for a more physiologically realistic occlusion, achieving a controlled compromise between contact stability and non-penetration (
Figure 6).
Ultimately, as the latter method resulted in the most consistent with the physiological needs of a stable occlusion, it has been adopted for the subsequent steps of this study. However, the framework has been designed to enable the practitioner to choose between the two approaches.
Manual stepwise refinement GUI
A GUI for potential final manual refinement of the occlusion obtained after the automatic phase was implemented.
To enhance the visualization of the dental arches, the complete models obtained from the IOS were transferred to the current occlusal position defined on the depth maps. A rigid registration procedure based on fiducial points was implemented. This procedure requires the operator to position at least three pairs of corresponding fiducial points to ensure reliable spatial alignment [
42] (
Figure 7).
Specifically, as shown in
Figure 7, the points were selected considering the cusps of the following teeth:
mesiobuccal cusps of the right and left second molar;
right and left canine cusps;
the midpoint of right and left maxillary central incisors’ incisal edge.
To improve local accuracy, this initial transformation is refined using traditional ICP together with a multi-resolution strategy that refines the alignment by subsampling the point clouds, and a point-to-plane metric that improves convergence behavior.
Subsequently, the same registration and alignment procedure was applied to the mandibular arch.
At the end of registration procedures, models in final occlusion and with complete geometries to be directly used in preoperative surgical planning were obtained (
Figure 8).
At this point, the maxillary and mandibular arches are displayed in an interactive interface with adjustable display parameters. The mandible is held fixed as a reference, while the maxilla is interactively modified through rigid transformations, including translations in mm and rotations in degrees along all axis, controlled via sliders. These controls are expressed using clinically meaningful terminology and indications to ensure intuitive interaction (
Figure 9).
The GUI includes tools to restore the configuration, undo the last transformation, view the applied transformation values, and save the result.
Real-time quantitative and qualitative feedback is provided during the interaction. Contacts are visually highlighted, while a heat map dynamically represents the depth of penetration.
The non-penetration constraint can be deactivated in cases where premature contacts prevent the achievement of a stable occlusion. In such situations, the areas of interpenetration are highlighted and identified as regions requiring preoperative occlusal adjustment.
Complete workflow validation
For ensuring a reproducible validation procedure, controlled misalignments were artificially introduced in the MrO to simulate malocclusion conditions and evaluate the robustness of the full digital workflow with respect to a known reference standard.
A surgeon from the Maxillofacial Surgery Unit at AOU Città della Salute e della Scienza di Torino executed the complete digital procedure. After completion of the automatic phase, the surgeon was asked to evaluate whether any further correction was required based on the resulting occlusal relationship.
In such cases, a strict step-by-step protocol was provided, specifying the sequence of actions to be performed in an orderly manner. All adjustments were carried out by moving the maxillary arch while maintaining the mandibular arch as the fixed reference.
Initially, the surgeon was instructed to position the models in frontal view and proceed through the revision steps reported in
Table 1. Each clinical parameter had to be evaluated and, when considered necessary, corrected before proceeding to the subsequent step.
If satisfied with the result obtained in the frontal view, the surgeon was instructed to proceed rotating the model in the right lateral view (
Table 2).
After all these ordered evaluations and corrections had been completed, and the surgeon was satisfied, resulting occlusal model was saved and stored as an .STL file.
In this validation phase, a transformation matrix was calculated to describe the rotational and translational discrepancies between the occlusion obtained using the full digital workflow and the MrO, with the latter being considered the reference standard. For the whole dataset, the matrix was calculated first after the automatic phase alone, and then after the manual stepwise refinement.
Statistical Analysis
Summary statistical indicators were calculated (Microsoft Excel, version 16.109; Microsoft Corp., Redmond, WA, USA) to facilitate the comparison of the obtained results:
Mean Value (MV), which represents the average value of the observed errors across all patients;
Mean Absolute Value (MAV), which represents the average magnitude of the deviations regardless of their sign;
Standard deviation (SD), to quantify the variability of the data from the mean value.
For each translational (X, Y, and Z axis) and rotational (pitch, roll, and yaw) discrepancy, values obtained after the automatic alignment phase were compared with those recorded after the manual refinement phase to determine whether the refinement procedure significantly reduced positioning errors. Results are reported as mean ± standard deviation (SD), together with the corresponding p - values.
Before comparative analysis, the normality of the paired differences was assessed using the Shapiro–Wilk test. Variables with normally distributed paired differences were analyzed using paired Student’s t-test, whereas Rotation on the Y - axes, which showed a non-normal distribution, was analyzed using the Wilcoxon signed-rank test. Statistical significance was set at p < 0.05.