Automatic Method for Bone Segmentation in Cone Beam Computed Tomography Data Set

Abstract

Due to technical aspects of Cone Beam Computed Tomography (CBCT), the automatic methods for bone segmentation are not widely used in the clinical practice of endodontics, orthodontics, oral and maxillofacial surgery. The aim of this study was to evaluate method’s accuracy for bone segmentation in CBCT data sets. The sliding three dimensional (3D) window, histogram filter and Otsu’s method were used to implement the automatic segmentation. The results of automatic segmentation were compared with the results of segmentation performed by an experienced oral and maxillofacial surgeon. Twenty patients and their forty CBCT data sets were used in this study (20 preoperative and 20 postoperative). Intraclass Correlation Coefficients (ICC) were calculated to prove the reliability of surgeon segmentations. ICC was 0.958 with 95% confidence interval [0.896 ... 0.983] in preoperative data sets and 0.931 with 95% confidence interval [0.836 ... 0.972] in postoperative data sets. Three basic metrics were used in order to evaluate the accuracy of the automatic method—Dice Similarity Coefficient (DSC), Root Mean Square (RMS), Average Distance Error (ADE) of surfaces mismatch and additional metric in order to evaluate computation time of segmentation was used. The mean value of preoperative DSC was 0.921, postoperative—0.911, the mean value of preoperative RMS was 0.559 mm, postoperative—0.647 mm, the ADE value of preoperative cases was 0.043 mm, postoperative—0.057 mm, the mean computational time to perform the segmentation was 46 s. The automatic method showed clinically acceptable accuracy results and thus can be used as a new tool for automatic bone segmentation in CBCT data. It can be applied in oral and maxillofacial surgery for performance of 3D Virtual Surgical Plan (VSP) or for postoperative follow-up.

1. Introduction

Three dimensional (3D) segmentation of bones from Cone Beam Computed Tomography (CBCT) is very important for the orthodontics, endodontics, oral and maxillofacial surgeons. Accuracy of segmentation ensures a correct diagnosis, an accurate 3D Virtual Surgical Plan (VSP) and a successful postoperative follow-up of the patients with Cranio-maxillofacial (CMF) deformities [1,2]. Application of CBCT imaging modality in oral and maxillofacial surgery is increased due to its lower radiation dose and shorter acquisition time compared with the conventional multislice CT (MSCT). However, images of CBCT most often are noisier and have beam hardening artifacts. Properties of CBCT also affect the ability to correctly display Hounsfield units (HU) of head tissues (immediately after surgery when edema of soft tissues is present) as opposed to traditional CT [3,4]. These disadvantages also influence image quality and accuracy of bone segmentation [5]. Due to these reasons, 3D VSP and evaluation of postoperative follow-up must be performed using segmentation by highly experienced experts. Such a procedure of bone segmentation becomes time-consuming [6].

Wang et al. published three studies dedicated to fully automatic bone segmentation using CBCT images [6,7,8]. In 2013 [6] and 2014 [7] studies, the principle of the automatic method was a patch-based sparse representation. Patient-specific atlas (probability map) using a sparse label fusion strategy from conventional CT atlases was used as a first estimation. Then a convex segmentation framework was used in order to get the final result. The introduced method led to an accurate segmentation but the basic limitation of this method was computation time (about 5 h [7]). Also, the variety of CBCT data sets (patients with metallic implants, metallic plates and different face deformities) was low in their studies. Data collection of CT to get more variety of the atlases is complicated in relation to the requirements of bioethics. In a 2016 study [8], the same authors also suggested a new automatic method that used the random forest. The multiclass classifier was used to create probability maps for each region of interest (mandible, maxilla and background). The results of the method were achieved almost the same as in the previous study [7] but also with some similar limitations: the limited amount of CBCT data and also computation time of segmentation (20 min). Gollmer and Buzug presented a fully automatic method for mandible segmentation. Segmentation was based on the idea of using a statistical shape model (SSM). They showed accurate results, however, they also—as did the previous authors—tested the algorithm with just six CBCT data sets [9]. Fan et al. [10] proposed an automatic method for mandible segmentation. Marker-based watershed transform method was used in their study. The authors performed accurate and fast enough (12-14 min per data set) segmentation on 20 CBCT data sets. The errors of segmentation were obtained mostly in these three basic areas—around the wisdom teeth, condyles and dental enamel. The reasons of segmentation errors were the different (bigger) or the same intensity of selected basic markers (mandible and background). Performance of manual editing in these areas was recommended. Eijnatten et al. [11] carried out a review of literature on different methods for bone segmentation. The authors discovered that global thresholding is the most used method for bone segmentation. However, a limitation of this method is that the manual post-processing is required.

This study presents a fully automatic 3D thresholding method which is based on local statistical information. In the previous pilot study [12], the automatic segmentation of ten CBCT data sets of preoperative patients was performed. We increased the amount and variety (postoperative cases) of data sets (n = 40) in the current study. Compared with the previous study, we included more anatomical areas of interest in order to better evaluate the accuracy of the presented method. Otsu’s method [13], histogram filter and 3D sliding window [14] are used for the implementation of the segmentation method. The accuracy of the current method is evaluated using three metrics—Dice Similarity Coefficient (DSC) [15], Root Mean Square (RMS) [16,17], Average Distance Error (ADE). The computation time of the segmentation is also evaluated. The results of the automatic segmentation are compared with the segmentation results of an experienced oral and maxillofacial surgeon.

2. Materials and Methods

2.1. Data Acquisition

A retrospective study was performed by using 40 CBCT data sets from the database of the Simonas Grybauskas’ Orthognathic Surgery clinic. Before the study, all CBCT data sets were anonymised in order to protect the patients’ data. Half of them (n = 20) were preoperative (obtained one week before the surgery) and another half (n = 20) were postoperative (obtained about one week after surgery) scans of the same patients’ group. All scans were done using the i-CAT FLX V17 machine. All patients undertook double jaws correction. CBCT data sets were acquired with a resolution of an isotropic voxel of 0.3 mm, 230 mm × 170 mm field of view (FOV), the time of exposure was 7 s, tube voltage 120 kV and tube current 5 mA. Preoperative and postoperative CBCT data sets were segmented by an experienced oral and maxillofacial surgeon using ITK-SNAP software (Version 3.4.0) selecting a global threshold value for each case individually. Mandible and lower parts of skull (including maxilla, zygomatic bone) were used to perform the segmentation. Therefore, mentioned anatomical regions were selected as the segmentation target assessing the importance of them to perform a surgery. The study framework is presented in Figure 1.

2.2. Description of Proposed Method

The proposed method was based on the finding of an optimal threshold by using Otsu’s method in a sliding 3D window. However, Otsu’s method works most accurately when the analyzing histogram is bimodal. This means that the intensity values of voxels must be divided into two basic classes $C_0$ with intensity range [1, …, T] and $C_1$ with intensity range [T + 1, …, L] (where L is the upper limit of intensity in the volume), T representing the threshold optimally separating modes in bimodal histogram. The number of voxels with intensity threshold i is denoted by $n_i$. A probability of intensity threshold i in an image is

$$p_i={\displaystyle \frac{n_i}{N}}$$

(1)

N—the total number of voxels. Then, the probabilities that a randomly selected voxels belong to one of the classes $C_0$ or $C_1$ are

$$C_0:{\omega }_0\left(T\right)=\sum _{i=0}^T{p}_i$$

(2)

$$C_1:{\omega }_1\left(T\right)=\sum _{i=T+1}^L{p}_i.$$

(3)

The average of classes is defined

$$C_0:{\mu }_0=\frac{{\displaystyle \sum _{i=1}^T}i·p_i}{{\omega }_0\left(T\right)}$$

(4)

$$C_1:{\mu }_1=\frac{{\displaystyle \sum _{i=T+1}^L}i·p_i}{{\omega }_1\left(T\right)}.$$

(5)

The intensity average of the total image is defined by

$${\omega }_0{\mu }_0+{\omega }_1{\mu }_1={\mu }_T$$

(6)

$${\omega }_0+{\omega }_1=1.$$

(7)

Using discriminant analysis, Otsu’s method defines both classes variance of the thresholded image

$${\sigma }_B^2={\omega }_0{({\mu }_0-{\mu }_T)}^2+{\omega }_1{({\mu }_1-{\mu }_T)}^2.$$

(8)

Then the optimal threshold value is calculated by

$$T^{*}=ArgMax\left\{{\sigma }_B^2\left(T\right)\right\}$$

(9)

$$1\le T<L.$$

(10)

The illustration of the proposed method is presented in Figure 2.

The typical histogram of CBCT data set is not bimodal and Otsu’s method does not work accurately in order to find optimal threshold value. The typical histogram of CBCT data set is presented in Figure 3.

In order to perform an accurate segmentation using Otsu’s method, a histogram filter was involved. The purpose of the filter is to remove irrelevant areas and to make analysed histograms bimodal. The optimal threshold was found after filtration, that is, removal of the $1\mathrm{st}$, $2\mathrm{nd}$ and $5\mathrm{th}$ areas in the histogram (Figure 3). Intensity values for histogram filtering areas were defined by the results of studies published in References [18,19]. Misch et al. [18] determined that the intensity of cortical bone of mandible is 1250 HU and it cover about 8% of the volume of mandible, trabecular bone intensity of mandible is 850–1250 HU, 350–850 HU levels corresponds to trabecular bone intensity of maxilla, levels of 150–350 HU limit posterior maxilla and levels of 150 HU could be found in sinus areas where the bone is very soft. The study performed by Norton et al. [19] divided the intensity scale of bone into these ranges of HU: anterior mandible 850 HU, posterior mandible/anterior maxilla 500–850 HU, posterior maxilla 0–500 HU, tuberosity region 0 HU. Air intensity is −1000 HU and internal anatomical areas such as cavity of maxillary sinuses, nasal cavities, airways are between −950 HU and −500 HU. This analysis was used to set which areas of the histogram are useful for the segmentation and which areas should be filtered out. Then the filter was built to automatically remove areas in the sliding volume histogram with the intensities of more than 1250 HU and areas with the intensities less than −500 HU. An optimal threshold value using Otsu’s method is found in the defined range [−500 HU ... 1250 HU]. The filter is involved before processing CBCT data set and it ensures that CBCT volumetric histogram is bimodal (Figure 4, C0—object, C1—background). Then the matrix of local thresholds (T(x,y,z)) is filled by the found T values and the final segmentation (A(x,y,z)) is performed by:

$$\mathbf{A}(\mathrm{x},\mathrm{y},\mathrm{z})=\left\{\begin{array}{c}0\mathrm{if}\mathbf{I}(\mathrm{x},\mathrm{y},\mathrm{z})<\mathbf{T}(\mathrm{x},\mathrm{y},\mathrm{z})\hfill \\ 1\mathrm{otherwise}.\hfill \end{array}\right.$$

(11)

Here A(x,y,z)—3D final matrix of segmentation, I(x,y,z)—3D matrix of the original CBCT data set.

The local threshold value T is found in a volume of sliding window (S(x,y,z)), the segmentation object is the segment of bone (mandible and maxilla)—C0 and background is the segment of soft tissues (including soft bone)—C1 (Figure 4).

2.3. Surface Reconstruction

Segmented bone (A(x,y,z)) is saved as a surface in the STereoLithography (STL) file format. The reconstruction of the bone surface from segmented voxels is performed by the volumetric reconstruction algorithm (Visualization and Computer Graphics (VCG) reconstruction filter) using MeshLab software (Visual Computing Lab of ISTI–CNR, University of Pisa) [20,21] (Figure 5).

The basic parameters of the surface reconstruction ((1) voxel size, (2) the level of the subvolume reconstruction process, (3) geodesic weighting, number of Volume Laplacian iterations, (4) widening, (5) number of smoothing iterations) were used the same for all forty cases in order to have comparable results of reconstruction.

2.4. Evaluation of Method Accuracy

The results of the proposed algorithm were compared with the results of segmentation performed by an experienced surgeon. In order to evaluate the reliability of surgeon’s segmentation, the segmentations were done twice with a 2 weeks interval. For the quantitative evaluation of the reliability, Intraclass Correlation Coefficient (ICC) by using a two-way mixed model, unit: single rater and the type of relationship: absolute agreement method [22] was calculated by:

$$ICC=\frac{M{S}_R-M{S}_E}{M{S}_R+(k-1)M{S}_E+\frac{k}{n}(M{S}_C-M{S}_E)},$$

(12)

where $M{S}_R$ = mean square for rows, k = number of raters/measurements; $M{S}_E$ = mean square for error; $M{S}_C$ = mean square for columns; n = number of subjects.

ICC was calculated separately for preoperative and postoperative cases and segmentation thresholds selected by the surgeon were used for this.

The evaluation of segmentation accuracy by using proposed method three basic metrics were used:

(1)

Root mean square (RMS) of the intersurface distance used to evaluate reconstructed surface mismatch:

$$RMS={\displaystyle \frac{1}{\sqrt{n}}}·\sqrt{\sum _{x,y,z=1}^n{(a_{x,y,z}-b_{x,y,z})}^2},$$

(13)

where $a_{x,y,z}$ represents the coordinates of reference surface point, $b_{x,y,z}$—the coordinates of surface point created with the proposed method;

(2)

Dice similarity coefficient (DSC), used to evaluate volume discrepancy:

$$DSC=\frac{2|A\cap B|}{\left(\right|A|+|B\left|\right)},$$

(14)

where A represents the volume of the reference model, B—volume of automatically segmented bone;

(3)

Average intersurface distance error (ADE) was calculated by:

$$ADE=\frac{1}{N}·\sum _{i=1}^n(a_{x,y,z}-b_{x,y,z}),$$

(15)

where $a_{x,y,z}$ represents the coordinates of 3D point of the reference model, $b_{x,y,z}$—the coordinates of 3D point of the automatically segmented model;

Positive and negative ADE values were calculated also.

Additional metric to evaluate the rapidity of segmentation the computation time of automatic segmentation was calculated. In this study, personal computer with parameters—processor—Intel(R) Core(TM) i7-4790 CPU @ 3.60 GHz, RAM—16 GB, system type—64-bit Windows 10 Operating System was used. The automatic method was implemented by using Matlab software.

3. Results

The reliability of segmentations performed by surgeon was evaluated by calculating ICC values for preoperative and postoperative CBCT data sets (

Table 1. Results of Intraclass Correlation Coefficient (ICC) using single rater, absolute agreement and 2-way random effects model.

	Itraclass Correlation	95% Confidence Interval		F Test with True Value 0
		Lower Bound	Upper Bound	Value	df1	df2	Sig
Single measures preoperative	0.958	0.896	0.983	49.03	19	19	0.000
Single measures postoperative	0.931	0.836	0.972	27.43	19	19	0.000

).

The results show that in preoperative data ICC = 0.958 with 95% confidence interval [0.896 ... 0.983] and in postoperative data ICC = 0.931 with 95% confidence interval [0.836 ... 0.972]. Calculated values of ICC show that the level of surgeon reliability is sufficient [22].

(1) The mean RMS value of intersurface distance in preoperative cases was 0.559 mm and in postoperative—0.647 mm. Calculated RMS values of intersurface distance for all cases are presented by a boxplot in Figure 6.

The interquartile range of boxplot is narrower for preoperative RMS values. The bigger distribution of RMS values is seen in postoperative cases. This could have been caused by a better quality of preoperative CBCT data sets.

(2) The mean value of DSC was 0.921 in preoperative cases. In postoperative cases, the mean value of DSC was 0.911. Calculated DSC values are similar and are very high—more than 0.9. The results of DSC in preoperative and postoperative segmentation data are presented by using the boxplot function in Figure 7.

The narrower range of interquartile is found in postoperative cases. Two DSC values are out of the DSC range in postoperative cases and are marked as outliers.

(3) The ADE values were calculated and divided into three groups: general average, positive average and negative average. The whole results are presented in

Table 2. The comparison of different automatic methods by 4 metrics.

Metric	Wang et al. [6]	Wang et al. [7]	Wang et al. [8]	Fan et al. [10]	Proposed
RMS, mm	-	-	-	-	$0.559{}^a$|$0.647{}^b$
DSC	$0.91{}^c$|$0.87{}^d$	$0.92{}^c$|$0.87{}^d$	$0.94{}^c$|$0.91{}^d$	0.97	$0.921{}^a$$0.911{}^b$
ADE, mm	0.61	0.65	0.42	-	$0.043{}^{1a}$|$0.279{}^{2a}$|$-0.242{}^{3a}$
					$0.057{}^{1b}$|$0.293{}^{2b}$|$-0.282{}^{3b}$
Time, s	-	18,000	1200	720–840	$45{}^a$|$47{}^b$

a—$preoperative$, b$-postoperative$, c—$region$$of$$interest$$maxilla$, c—$region$$of$$interest$$mandible$, 1—$average$$distance$$error$, 2—$possitive$$average$$distance$$error$, 3—$negative$$average$$distance$$error$.

.

The mean value of computational time to perform an automatic segmentation by using all CBCT data sets was 46 s. Compared with other studies, this computation time is very fast. The achieved results are compared with the results of other similar studies. The comparison is presented in Table 2.

DSC is the most popular metric in order to prove the accuracy of segmentation. Achieved values of DSC were similar in comparison with other studies (values were near 0.9). ADE also is a frequent metric in order to evaluate the segmentation. In our study, ADE was calculated and divided into three groups. The results between ADE values in different groups (preoperative and postoperative) were compared. It was found that differences were similar (Table 2). In comparison with the other studies, in which the proposed method was used, ADE values were found lower. The proposed method showed a short computation time to an automatic segmentation performance. The comparison of segmentation results by using automatic and surgeon methods are shown in Figure 8.

The results of segmentation by using the proposed method showed that bones with low density were not fully segmented. However, areas of condyles and sinuses were segmented with fewer holes by using the proposed method, compared to the surgeon segmentation results (Figure 8).

4. Discussion

Automatic methods for bone segmentation are very important in order to get 3D surfaces. It helps surgeons in making correct diagnoses, in performing accurate and quicker VSP and in the evaluation of postoperative follow-up without the influence of surgeon’s experience [23,24,25,26,27]. The aim of this study was to evaluate method’s accuracy for bone segmentation in CBCT data sets. The method is based on the histogram filter, 3D sliding window and Otsu’s thresholding. The results of automatic segmentation revealed sufficient accuracy of bone segmentation. For the evaluation of method accuracy, three kind of metrics were used: voxel-based–DSC [15,17,28], based on the intersurface distance evaluation–RMS and ADE [16,17]. Additional metric based on the time to perform segmentation–computation time was calculated. The mean DSC coefficient values of two groups (preoperative and postoperative) were bigger than 0.9, which proved the complete overlap and high accuracy of the automatic and surgeon segmented bone voxels. The mean RMS values of intersurface distance for preoperative (0.559 mm) and postoperative (0.647 mm) cases were about two times bigger than the voxel size (0.3 mm). Calculated ADE values showed small discrepancies of surfaces. It shows that segmentation result is good, either we succeeded to avoid superimposition step because both automatic and global segmentations were made by using the same source data sets. In this way, surface superimposition did not yield any additional errors. Compared with the other studies [7,8,10], the proposed method performed the segmentation very rapidly ( 46 s/case). The achieved results showed that the proposed automatic method worked accurately.

However, the proposed method had some limitations. CBCT images were not filtered for metal artefacts (brackets, metal plates and mini-implants) before or after segmentation. Due to this reason, metal artefacts were seen in 3D models. These artefacts hide important areas of bone. Therefore, assessment of the bone near these artefacts became complicated and inaccurate. Especially it is important when postoperative follow-ups are performed [29,30]. The other limitation of the proposed method is the difficulty to segment areas with low density (thin anatomical areas, e.g., alveolar part of the mandible, mandibular condyles or areas of maxillary sinuses) in CBCT images [31]. This is also valid for other threshold-based segmentation methods. Problematic areas could be segmented including more sensitive techniques [31,32,33,34].

Further studies may be concentrated to evaluate proposed method with higher amount and different kind of CT/CBCT data sets. The next direction of further study may be to increase accuracy of the segmentation for anatomical areas with low density. Fully automatic segmentation of selected anatomical areas especially condyles would be important tool to increase the evaluation of treatment or postoperative follow-up [35,36].

5. Conclusions

The presented study proposed a new automatic method for bone segmentation in a CBCT data set. The important feature of the proposed segmentation method is simple and rapid implementation. The results of segmentation were very accurate and reliable to use in the clinical practice of oral and maxillofacial surgery. The method does not require access to a computer with high computation power parameters. It can be integrated into the most popular software of 3D medical imaging processing using an ordinary computer.