1 INTRODUCTION

Unilateral cleft lip and palate (UCP) is a common congenital maxillofacial hypoplasia that exhibits a high incidence. Data for 7.5 million births from the international perinatal database of typical oral clefts (IPDTOC) showed that the prevalence of cleft lip and palate was 6.64 per 10 000 births worldwide.1 The resultant deformity, located in the maxilla and other midfacial areas, is characterized by the incomplete formation of the lip, alveolar crest, hard palate and soft palate.2 It often results in some combination of feeding, deglutition, speaking, hearing and/or psychological problems.3

Facial morphology plays a crucial role in mental health and social recognition. The fundamental goal of reconstruction is to restore the facial symmetry three-dimensionally (3D) to within a clinically acceptable range of the general population.4 Thus, there is an increasing need to quantitatively measure and accurately describe the extent of facial asymmetry. Many studies have attempted to investigate the morphological features of the maxilla in UCP patients.4-7 There is still no consistent conclusion in terms of whether significant asymmetry exists in the “deeper” maxilla. Several studies demonstrated that the maxillary asymmetry is confined to “local” regions of the cleft,5, 6 while others stated that the asymmetry also involves “deeper” regions of the maxilla.4, 7 From an etiological point of view, no single factor has been identified as a cause of the maxillary asymmetry. Several studies indicated that numerous etiological factors might contribute to this asymmetry, such as the location and extent of the defect, deviation of the nasal septum cartilage and side effects of previous surgery.8, 9 Since the cleft is located in the maxilla, it is imperative to understand the morphology of the defect and the defect-related maxillary dysmorphology. Barbosa et al10 revealed some clues indicating that the defect volume is related to the gap, arch, nasal and dental parameters to a certain degree. However, few studies focused on revealing the morphological relationship between the defect and affected maxillary half.

An essential step in the morphological analysis workflow is to obtain a virtual 3D model of the maxilla and defect. Although manual segmentation is considered the gold standard, it has the inherent drawback of being time-consuming, which affects its application in large-scale clinical studies. With the rapid advance in artificial intelligence (AI), several methods for auto-segmentation of the region of interest have been proposed in the literature.11, 12 Chang et al13 proposed a model-based segmentation approach for a particular interest in the outer surface of the anterior wall of the maxilla. However, the thin bony structure surrounding the maxillary sinus is challenging to identify and segment from the adjacent maxillofacial bones. Compared with the auto-segmentation of the maxilla, the segmentation of the defect is more difficult. The defect has no specific boundary with the nasal and oral cavities, and its image intensity is similar to the soft tissue in the CBCT. Recently, Zhang et al16 utilized 3D U-Net14, 15 which incorporated non-rigid volumetric registration to explore the discrepancy between complete and diseased maxillae. Using this algorithm, they obtained more accurate segmentations of the maxilla and cleft for UCP patients, compared to other methods. Despite the fact that these state-of-the-art AI algorithms still cannot be directly applied in a clinical study in a fully automatic manner, we can use them to obtain initial segmentation results, followed by manual post-processing refinement, to achieve clinically applicable segmentations on our dataset. In this way, the preparation time for a large sample size can be significantly reduced.

This study aimed to investigate the defect factors responsible for the variability of the maxilla on the cleft side by quantifying the 3D asymmetry of the maxilla in patients with UCP based on a deep learning-based CBCT image auto-segmentation method with manual refinement.

2 MATERIALS AND METHODS 2.1 Subjects

This retrospective study was performed under the approval of the institutional review board of Beijing Tian Tan Hospital, Capital Medical University (KY2017-072-01). The CBCT images of 60 subjects were selected from the Department of Orthodontics, Beijing Stomatological Hospital, Capital Medical University. All of the subjects had been diagnosed with non-syndromic UCP and received primary lip and palate repair. The subjects who had a history of alveolar bone grafts, orthodontic treatment, maxillofacial neoplasia, trauma or orthognathic surgery were excluded. A total of 60 subjects were included in this study, consisting of 39 males and 21 females, with a mean age of 11.52 years (SD = 3.27 years; range of 8-18 years), presenting with 41 left-side defects and 19 right-side defects.

All of these CBCT images were obtained during routine radiographic documentation prior to orthodontic treatment or surgical planning. All CBCT images were acquired on the NewTom VG scanner (QR srl, Verona, Italy) under a standard scanning protocol: 110 kV, 6.35 mA, 15 × 15 cm field of view and 0.250 mm slice thickness. The subjects were in an upright position for maximum intercuspation. The Frankfort plane was parallel to the floor.

2.2 Maxilla and defect segmentation

The 60 CBCT image sets were exported in DICOM format. To reduce the learning difficulty for the AI algorithm, the orientation of each CBCT image was adjusted to have every defect on the left side of the maxilla. We randomly selected 30 CBCT images and included them in Group 1. The remaining 30 images were assigned to Group 2. All the CBCT images in Group 1 were annotated manually and semi-automatically using ITK-SNAP17 (version 3.6.0; www.itksnap.org) by two trained specialists, serving as the gold standard. The manual segmentations from the two specialists had Cohen's kappa statistic of 0.98, indicating our manual segmentations were in the high standard. For images in Group 1, the maxilla was manually segmented by tracing the boundaries where it articulated with the adjacent bony structures in the craniofacial region. The defect was manually segmented by following the contralateral shape to obtain a continuation of the alveolar ridge and hard palate.

The automatic segmentation protocol used in this study is 3D U-Net,14, 15 which is a well-known deep learning neural network for biomedical image segmentation. The network architecture of 3D U-Net consists of three elements: contracting path, expansive path and skip connections between the two paths. The contracting path is used to capture the context in the image, whereas the expansive path is used to recover local information based on deconvolution. The skip connections between the two paths help to fuse the local and global features, resulting in its outstanding performance in image segmentation.

A 3D U-Net model was trained and tested based on the 30 CBCT images in Group 1, where 24 images served as the training set, three images as the validation set, and the remaining three images as the test set. In the training set, each CBCT image generated 3,000 patches with a voxel size of 64 × 64 × 64 with voxel spacing of 0.25 mm, which served as the real training samples. The amount of training data in terms of patches was much larger than the number of original images, which helped to overcome the limited number of samples in the biomedical field. The model was optimized to have minimal generalized Dice loss18 for 50 epochs using the Adam optimizer.19 All processes were implemented using PyTorch,20 an open-source deep learning library developed by Facebook.

Each CBCT was split into several patches with the same size of training patches in the sliding window manner. The output was the patch containing the voxel label, and its size was the same as the input patch. Combining the output patches, we obtained the segmentation result of the CBCT. The accuracy of the segmentation model is evaluated by the Dice similarity coefficient (DSC), defined as follows: where and represent the cardinalities of the learned and manual sets, and represents the intersection of the two sets. A value of 0 indicates no similarity, whereas a value of 1 indicates perfect agreement.

The trained 3D U-Net model was then used to predict the segmentation for the remaining 30 CBCT images in Group 2. An orthodontist further refined the predicted segmentation results following the automatic segmentation procedure. In addition, to eliminate the effect of the number of teeth on the volumetric measurements (especially for cases at different stages of dental development), all crowns were manually removed in all of the CBCT images (both Group 1 and Group 2). All 60 pairs of segmented maxilla and defect models were then used for morphometric and statistical analysis.

2.3 Description of measurement

For 3D morphometric quantification, four midsagittal landmarks and 10 bilateral landmarks, as defined in Table 1 and illustrated in Figure 1A,B, were identified on the surface of the 3D segmented model and verified in the multiple planar reformat mode. Three reference planes (the Frankfort horizontal plane, midsagittal plane and coronal plane) established a coordinate system (Table 1). The maxilla was separated by the midsagittal plane for measurements on the cleft and non-cleft sides. Eleven structural parameters of the maxilla (10 linear distances and 1 volume, including the maxilla and alveolar crest) were measured in horizontal, midsagittal and coronal plane projections to characterize the maxillary asymmetry. In addition, four structural parameters of the defect (three linear distances and one volume) were measured, in order to investigate their correlation with those of the maxilla on the cleft side. All parameters are defined in Table 1, while various defect measurements are shown in Figure 1C.

TABLE 1. 3D cephalometric landmarks, reference planes and measurements of the maxilla and defect Items Definition Landmarks N Intersection of internasal suture with nasofrontal suture S Midpoint of sella turcica ANS Most anterior point of anterior nasal spine PNS Most posterior point of posterior nasal spine Po Uppermost point on bony external auditory meatus Sm Superior most extent of the maxilla Or Lowest point on infraorbital edge Lap Most lateral points on the nasal aperture J Intersection of the outline of the tuberosity of the maxilla and zygomatic buttress Mt Posterior most extent of the maxillary tuberosity Am Posterior most extent of the anterior contour of the maxilla Spc Midpoint of labial alveolar crest of maxillary canine (No missing canine was observed in all subjects) Spm Midpoint of buccal alveolar crest of maxillary first molar Aa Most inferior anterior point of alveolar crest Reference planes Horizontal plane (FH plane) Plane that passes through the bilateral Po and Or on the non-defect side Midsagittal plane (MS plane) Plane perpendicular to the FH plane passing through the N and S Coronal plane (CR plane) Plane perpendicular to the FH and MS plane passing through the N Measurements Maxillay length (Lmax) Sagittal distance from Am to Mt Maxillary anterior width (AntWmax) Transverse distance from Lap to the MS plane Maxillary posterior width (PosWmax) Transverse distance from J to the MS plane Maxillary anterior height (AntHmax) Vertical distance from Or to ANS Maxillary posterior height (PosHmax) Vertical distance from Sm to PNS Maxillary volume (Vmax) Volume of the segmented individual maxilla Alveolar length (Lalv) Maximum sagittal distance from Aa to Mt Alveolar anterior width (AntWalv) Transverse distance from Spc to the MS plane Alveolar posterior width (PosWalv) Transverse distance from Spm to the MS plane Alveolar anterior height (AntHalv) Vertical distance from Spc to ANS Alveolar posterior height (PosHalv) Vertical distance from Spm to PNS Defect length (Ldef) Maximum Sagittal distance of the defect Defect width (Wdef) Maximum transverse distance of the defect Defect height (Hdef) Maximum vertical distance of the defect Defect volume (Vdef) Volume of the segmented defect

The main landmarks on the 3D segmented model. (A) Frontal view (cleft side and non-cleft side; landmarks on the non-cleft side are denoted by superscript '); (B) Lateral view (cleft side); (C) 3D segmentation measurements of the defect (Ldef, Wdef, and Hdef indicate the defect length, width, and height, respectively)

The lengths, widths and heights were measured by calculating the distances between the position (voxel coordinates) of landmarks, and the volumes of the maxilla and defect were measured based on the segmentation voxel counting. These measurements were carried out using ITK-SNAP.

2.4 Statistical analysis

To assess intra-observer consistency and inter-observer reliability, landmark identifications and distance measurements were repeated two weeks after the first examination by two trained specialists. The intra-observer consistency, which is measured by the intra-class correlation coefficient (ICC), is greater than 0.90, confirming the consistency of the measurements. The inter-observer reliability, measured by the difference between the two observers, was tested by a paired t-test, having the P-value of .39, indicating a good agreement between two observers. The mean values of measurements were used for statistical analysis. The data are presented as mean values and standard deviations. The maxillary asymmetry between the cleft and non-cleft sides was compared using the paired samples t-test. A multiple linear regression was carried out to analyse the relationship between the parameters of the defect and those of the cleft side of the maxilla, with adjustments for the age and gender of subjects in the regression. The Bonferroni correction was used to control for the family-wise type-I error probability of multiple comparisons. All statistical analyses were performed using SPSS (Version 19.0; IBM Co.). The level of significance was set at P < .05.

3 RESULTS 3.1 Maxilla and defect segmentation

The Dice similarity coefficients of the maxilla and defect between the manual and auto-segmentation samples were 0.92 ± 0.01 and 0.77 ± 0.06, respectively. Figure 2 shows the automatic (the first row) and manual (the second row) segmentations of three test samples in Group 1. The processing time of the automatic segmentation on one image was approximately one min/CBCT image set, using a GPU with a model of NVIDIA GTX 2080 Ti. Although the automatic results still needed an orthodontist to refine the images (average of approximately 5 min/CBCT image), the total processing time to obtain the final accurate segmentation results for Group 2 was significantly reduced, compared to the time for Group 1 (average of approximately 10 hours per CBCT images).

The segmentation results of the three test samples in Group 1. The first and second rows represent the results from 3D U-Net and manual, respectively. The red and orange parts represent the maxillae, and the green and yellow parts represent the cleft defects

3.2 Maxillary asymmetry analysis and defect measurements

Statistically significant differences were observed upon maxillary asymmetry analysis, shown in Table 2. An overall difference in the measurements was noted on the cleft side of the maxilla, mostly concerning the pyriform aperture and alveolar crest area. The cleft side of the maxilla demonstrated significantly reduced values of maxillary volume (Vmax) and length (Lmax) as well as alveolar length (Lalv), anterior width (AntWalv), posterior width (PosWalv), anterior height (AntHalv) and posterior height (PosHalv). A significant increase in maxillary anterior width (AntWmax) was observed for the cleft side of the maxilla. The defect volume (Vdef), length (Ldef), width (Wdef) and height (Hdef) had mean measurements of 1.24 ± 0.29 × 103 mm3, 18.76 ± 7.25 mm, 13.89 ± 1.89 mm and 13.40 ± 3.10 mm, respectively. The defect structure parameters are shown in Table 3.

TABLE 2. Measurement and analysis of the cleft and non-cleft sides of the maxilla Parameter Defect side Non-defect side P-value Mean SD Mean SD Vmax (×103 mm3) 18.02 3.24 19.32 3.53 .000* Lmax (mm) 34.31 2.56 41.27 3.72 .000* AntWmax (mm) 14.42 2.55 13.09 1.95 .002* PosWmax (mm) 34.73 2.91 35.08 2.44 .293 AntHmax (mm) 20.68 3.98 21.33 4.67 .236 PosHmax (mm) 32.91 4.31 32.89 4.33 .893 Lalv (mm) 36.81 3.95 42.27 5.25 .000* AntWalv (mm) 15.91 2.21 17.87 1.99 .000* PosWalv (mm) 27.80 2.10 28.79 1.83 .000* AntHalv (mm) 14.45 3.15 15.45 2.76 .000* PosHalv (mm) 10.39 3.19 11.13 3.38 .009* TABLE 3. Defect structure parameters and measurements, where Vdef denotes the volume of the segmented defect; and Ldef, Wdef, and Hdef denote the maximum sagittal, transverse and vertical distances of the defect, respectively Parameter Mean SD Vdef (×103 mm3) 1.24 0.29 Ldef (mm) 18.76 7.25 Wdef (mm) 13.89 1.89 Hdef (mm) 13.40 3.10 3.3 Multiple linear regression

Since cleft lip and palate can be caused by genetics, we considered defect parameters as independent variables and maxillary parameters on the cleft side as dependent variables. Each analysis was adjusted for age and gender. After performing a multiple linear regression followed by the Bonferroni correction, it was found that (a) the cleft side maxillary volume (Vmax), maxillary length (Lmax), alveolar anterior height (AntHalv) and alveolar posterior height (PosHalv) were significantly related to the defect height (Hdef); (b) the maxillary anterior width (AntWmax), alveolar anterior width (AntWalv) and alveolar posterior width (PosWalv) were significantly related to the defect width (Wdef); (c) the alveolar length (Lalv) was significantly related to the defect height (Hdef); and (d) the alveolar anterior width (AntWalv) was also significantly related to the defect volume (Vdef). The linear regression analysis results between the maxillary and defect parameters on the cleft side are presented in Table 4.

TABLE 4. Results of multiple linear regression analysis regarding the defect and relationship to the cleft side of the maxilla with adjusted age and gender Dependent variable Independent variable Coefficient Standard error P-value R2 adjusted Vmax (×103 mm3) Vdef (×103 mm3) 1.440 1.240 .252 .474 Ldef (mm) −0.067 0.044 .135 Wdef (mm) 0.307 0.190 .113 Hdef (mm) 0.329 0.108 .003* Lmax (mm) Vdef (×103 mm3) −1.304 1.156 .264 .272 Ldef (mm) −0.029 0.041 .485 Wdef (mm) 0.273 0.177 .128 Hdef (mm) 0.314 0.100 .003* AntWmax (mm) Vdef (×103 mm3) −1.487 1.224 .230 .180 Ldef (mm) −0.037 0.043 .394 Wdef (mm) 0.584 0.187 .003* Hdef (mm) −0.034 0.106 .747 PosWmax (mm) Vdef (×103 mm3) −1.652 1.375 .235 .201 Ldef (mm) −0.015 0.049 .757 Wdef (mm) 0.553 0.210 .011** Hdef (mm) 0.214 0.199 .077 AntHmax (mm) Vdef (×103 mm3) 1.263 1.556 .421 .512 Ldef (mm) −0.018 0.055 .749 Wdef (mm) 0.322 0.238 .182 Hdef (mm) −0.096 0.134 .480 PosHmax (mm) Vdef (×103 mm3) 3.059 1.957 .124 .265 Ldef (mm) −0.091 0.069 .194 Wdef (mm) −0.167 0.299 .580 Hdef (mm) 0.124 0.169 .466 Lalv (mm) Vdef (×103 mm3) −1.286 1.673 .446 .359 Ldef (mm) −0.132 0.059 .031** Wdef (mm) 0.526 0.256 .045** Hdef (mm) 0.427 0.145 .005* AntWalv (mm) Vdef (×103 mm3) −2.916 0.961 .004* .325 Ldef (mm) −0.044 0.034 .200 Wdef (mm) 0.615 0.147 .000* Hdef (mm) −0.053 0.083 .524 PosWalv (mm) Vdef (×103 mm3) 0.446 1.016 .662 .164 Ldef (mm) −0.023 0.036 .520 Wdef (mm) 0.477 0.155