A Fast Automatic Reconstruction Method for Panoramic Images Based on Cone Beam Computed Tomography

Abstract

Panoramic images have been widely used in the diagnosis of dental diseases. In the process of panoramic image reconstruction, the position of the dental arch curve usually affects the quality of display content, especially the completion level of the panoramic image. In addition, the metal implants in the patient’s mouth often lead the contrast of the panoramic image to decrease. This paper describes a method to automatically synthesize panoramic images from dental cone beam computed tomography (CBCT) data. The proposed method has two essential features: the first feature is that the method can detect the dental arch curve through axial maximum intensity projection images over different ranges, and the second feature is that our method is able to adjust the intensity distribution of the implant in critical areas, to reduce the impact of the implant on the contrast of the panoramic image. The proposed method was tested on 50 CBCT datasets; the panoramic images generated by this method were compared with images attained from three other commonly used approaches and then subjectively scored by three experienced dentists. In the comprehensive image contrast score, the method in this paper has the highest score of 11.16 ± 2.64 points. The results show that the panoramic images generated by this method have better image contrast.

1. Introduction

In the field of modern dentistry, cone beam computed tomography (CBCT) is currently a diagnostic imaging technique widely used in assisted diagnosis, virtual simulation, and treatment planning [1,2], with multiple clinical applications in oral implants [3,4,5,6]. The panoramic image displays tissue information, such as the maxillary, mandible, and dentition in one image by projecting the oral CBCT data within a certain range, providing the dentist with an intuitive image of the oral and maxillofacial tissues [7]; it plays an indispensable role in the post-processing of oral images. In addition, research on tooth segmentation based on panoramic images [8], mandibular neural canal segmentation [9], mandibular fracture detection [10], and oral 3D reconstruction [11] has gradually increased. Therefore, the automatic panoramic image reconstruction method can provide doctors with accurate images and provide a data basis for follow-up research.

In general, there are two reconstruction methods for oral panoramic images: the first method is to construct the image based on the dental arch line, through the curved planar reformation (CPR) [12], which directly generates the corresponding multi-plane reconstruction (MPR) image of the curve. This method can quickly reconstruct the image but may lead to the inability to fully display the root of the tooth, abnormal dentition, or implant in the upper and lower jaws [13,14]. The second method is to form the image based on the curve and thickness of the dental arch. The data within the thickness range of the dental arch are rearranged to produce the curved multiple MPR image sets, and the final oral panorama image is reconstructed by ray-summation or X-ray method [15,16], which displays the internal structure of the oral cavity on a higher completion level. However, this approach often brings unnecessary areas into the calculation process and raises the calculation cost, which results in reduced contrast and blurry images. Therefore, finding a small and suitable arch curve is the key to the reconstruction process.

In order to obtain accurate, high-quality panoramic images [17], related research has proposed a variety of methods for fitting dental arches. Chanwimaluang et al. [18] used a polynomial equation to fit the dental arch to generate a panoramic image from the dental arch curve. After that, more and more studies used different curves and control points [19,20] to fit the dental arch more accurately. To determine the shape of the dental arch in different patients, Bae et al. [21] used tooth position-controlled cubic B-splines to successfully fit different patient arch shapes. In terms of automatic reconstruction of panoramic images [22,23,24], Sa-Ing V et al. [25] introduced a dental arch detection algorithm based on threshold and image morphology algorithms. The dental arch calculation is simplified by selecting a slice of the original CBCT for downsampling. After the dental arch is detected, a series of curved MPR image sets are generated to produce a ray-sum panoramic image. In order to show more of the inside of the mouth in the panoramic image, Luo T et al. [26] described a method for reconstructing panoramic images based on three-dimensional surfaces. First, the dental arch curve was generated in axial maximum intensity projection (MIP), and the long axis curve was generated according to the shape of the long axis of the upper and lower teeth, extracting a three-dimensional surface through the curve to describe the entire dentition, and then expanding the three-dimensional surface to obtain a panoramic image. The panoramic images generated by this method can make the slices very thin to obtain sharper images and easily recognizable anatomical structures. In order to automatically obtain a high-contrast panoramic image, Yun Z et al. [27] introduced an automatic extraction method for dental arch curves. The coronal MIP histogram is fitted by the normal curve, the range of the transverse MIP section is determined according to the peak to detect the arch line and the arch thickness area, the curved MPR image set is generated in the arch thickness range, and finally, the panoramic image is reconstructed using the synthesis enhancement algorithm. All of these methods can automate the processing of oral CBCT data to reconstruct the panoramic image, but there are still problems with the contrast caused by metal implants in the panoramic image.

Based on previous research, we found that in the process of generating the dental arch curve, the dental arch can be divided into multiple regions according to the growth of the teeth, and different dental arches can be generated, respectively. To avoid shifting caused by anomalies and calculation errors, we propose a dental arch curve detection and adjustment method based on axial MIP images over different ranges. This method displays the entire dentition with a high level of completion in the panoramic image. On the other hand, many studies proposed different reconstruction methods to improve the contrast of panoramic images, while ignoring the contrast problems caused by metal implants in panoramic images. To resolve contrast issues, the strength of the metal implant is altered during the image extraction and fusion process to reduce the impact on other tissue areas in our method. By applying this technique, the contrast and quality of panoramic images were significantly improved.

The rest of this paper proceeds as follows. The full flow of the panoramic image algorithm is introduced in Section 2. In Section 3, the comparative experiments and evaluation results of different methods are presented, and the advantages and disadvantages of the proposed method are discussed. Finally, we make a conclusion in Section 4 and propose future directions. This method improves the contrast of panoramic images by altering the strength of the metal implant to reduce the impact on other tissue areas.

2. Materials and Methods

2.1. Data Acquisitions

The 50 clinical dental CBCT samples employed in this paper were collected from randomly selected patients by the Shanghai Fifth People’s Hospital and Zhongshan Hospital, Fudan University, Shanghai, China, and the involved equipment was acquired using a Planmeca ProMax 3D Max (Planmeca, Finland) scanner and NewTom VGi (QR s. r. l, Verona, Italy) scanner with the following parameters: 110 kVp, 3–8 mA, scan time: 27 s. The size of all slices of CBCT was 776 × 776 pixels, with a 0.2 × 0.2 mm pixel size. The inter-slice distance was 0.2 mm. The number of slices was around 432, depending on the patient situation.

2.2. Methods

As shown in Figure 1, the workflow of our proposed method was mainly divided into two steps. First, the axial MIP image of the mandible and the axial MIP image without the crown were calculated according to the coronal MIP, and the dental arch curve was detected; then, the enhanced image was extracted from the curved MPR image set, and the final panoramic image was obtained using the image fusion algorithm.

2.2.1. Dental Arch Curve Detection

At present, the main dental arch detection method is to detect the dental arch regarding the axial MIP. If the axial MIP is directly reconstructed using CBCT data, excessive amounts of irrelevant tissue will be displayed in the axial MIP (Figure 2a), which will result in poor accuracy of the dental arch detection. If a certain range of data is needed to reconstruct the axial MIP, the range is normally preset and manually modified by the user, which causes undesired interactions that cannot be automated. The MIP images generated from four different ranges of data are shown in Figure 2.

In the process of generating MIP images, the idea is to encode the volume data of the largest signal on the trajectory projected, along with the perspective of the coronal surface, and project the high-intensity signal of the three-dimensional space into a plane to form a continuous image.

In the coronary MIP image, the data on enamel, cortical bone, and soft tissue were the most numerous, and the Hu values of these three tissues were quite different, so they were represented as three independent peaks in the Hu intensity histogram, of which the largest peak at the abscissa was the distribution of the Hu value of the tooth. The tissue Hu value in the histogram can be obtained by fitting Formula (1). First, the generated Hu intensity histogram is smoothed using Gaussian smoothing to remove the violent fluctuations in the local areas of the histogram and reduce the impact on the crests. Then, the histogram curve is derived to satisfy Formula (2), where the first derivative is 0 and the second derivative is negative to determine the peak value and N. Finally, the final curve is obtained by fitting a Gaussian multimodal function [28].

$$fx={\displaystyle \sum }_{i=1}^N\frac{A_i}{\sqrt{2\pi {\delta }_i^2}}{e}^{-\frac{{\left(x-{\mu }_i\right)}^2}{2{\delta }_i^2}}$$

(1)

$$\frac{dp}{dx}=0, \frac{d^{2}p}{d{x}^2}<0$$

(2)

where $A_i$ is the peak gain corresponding to each peak, and the sum is equal to one; ${\mu }_i$ is the mean value corresponding to each peak; ${\delta }_i$ is the standard deviation corresponding to each peak. $N$ is the number of peaks of the Gaussian function; $p$ is the histogram curve.

After fitting the histogram with a Gaussian multimodal function, the threshold of the filtered image was obtained at the peak with the largest abscissa (Figure 3b). It can be seen from the coronal MIP image that the X-ray absorption of tooth enamel is higher than that of the upper and lower jaw cortical bones, and the intersection of teeth and alveolar bones can be more clearly distinguished. The distribution of the teeth can be obtained after a horizontal projection of the tooth by threshold filtering (Figure 3c), with the largest peak at the patient’s tooth bite near the upper and lower jaw boundaries. The slicing range values of crests and troughs, both oral occlusal and intersections of teeth and alveolar bone, are determined by the same process described above. Then, the axial MIP images are generated from two ranges of slices. This procedure of splitting the original image into two fragments for better imaging is named sub-regional detection, as shown in Figure 2c,d.

The calculation of the dental arch curve is established on the axial MIP. Bone thresholding parameters are calculated using the same method and provide a filtered mask image. In the mask image, the contour formed by the jaw and teeth is the largest in area, and small holes are commonly observed in the contour. Furthermore, visible burrs are often seen on the boundaries of contour, which leads to rough edges. Therefore, the contour pattern is treated by the expansion corrosion method to obtain a relatively full and smooth contour (Figure 4a). The set of contour points of all adjacent areas is extracted in the smooth contour map. Eventually, the contour with the largest area in the image is extracted, and the contour image of the mandible and teeth is obtained.

Detection of the dental arch curve based on mask images. Firstly, the mask image (Figure 4a) is processed by the thinning algorithm [29,30,31], and the bifurcation endpoint template of 8 fields is added to compare the skeleton line in the refinement calculation process (Figure 4b), the bifurcation and endpoint structure in the skeleton are recorded, and the branching statistics and fork removal are performed (Figure 4c). Then, the insertion point is determined on the skeleton. A recent study [32] shows that a smooth dental arch can be fitted by 11 insertion points with sufficient fitting accuracy and low computational cost. Through conducted experiments, 13 points were used to satisfy the fitting accuracy requirement. The insertion points are determined at equal intervals along the horizontal axis, generating a non-uniform cubic B-spline curve to represent the dental arch (Figure 4d). Finally, in order to make the insertion point fit the dental arch more accurately, we proposed an interpolation point calibration method. As shown by the yellow lines and red dots in Figure 5, at the insertion point within the range of teeth, a calibration line is made with the occlusal slice value, and the position of the insertion point is adjusted by calibrating the tooth center position on the line. The method can reduce the lack of fitting accuracy caused by the abnormal contour and the error of the refinement algorithm and can display the structure of the tooth and the implant more completely.

As can be seen from Figure 5, teeth in different positions have different degrees of inclination in space. For example, the molar area (Figure 5c) is nearly vertical; the crown in the cuspid area (Figure 5b) is slightly tilted forward, and the root is slightly backward; the teeth are significantly tilted in the incisor area (Figure 5c). If the curved MPR image is reconstructed based on a single arch curve, the internal condition of the tooth in the incision area is difficult to visualize, and if the curved MPR image is reconstructed based on the arch line and arch thickness, it is easy to add other non-target areas to participate in the reconstruction. We propose to reconstruct a curved MPR image set based on multiple arch curves. The distance of the dental arch curve transformation was calculated based on the skeletal lines (Figure 4c,f,g) of the axial MIP image of the mandible (Figure 2b) and the uncrowned axial MIP image (Figure 2c). The size of the curved MPR image set is determined by the length of the dental arch curve and the thickness of the CBCT data, respectively, so the curve needs to be isometrically transformed to generate a curved MPR image of consistent size and fused. The dental arch curve of the new position (Figure 4h) is obtained after the isometric transformation. Finally, using multiple arch curves to generate a curved MPR can show more roots or implants without adding too much computational cost.

2.2.2. Image Enhancement and Synthesis Algorithms

Enhanced images were generated from curved MPR image sets. In order to suppress non-interested tissues and regions, the parts of interest in curved MPR at different angles were extracted from each panoramic image, and their intensity distribution was amplified. We propose a nonlinear extraction algorithm. The enhanced image $P_0$ can be generated by the following formula:

$$P_0=\left\{\begin{array}{cc}\hfill \alpha Y\ast ln\left({{\displaystyle \sum }}_{n=1}^N{e}^{\raisebox{1ex}{$P_n\left(i,j\right)$}\!\left/ \!\raisebox{-1ex}{$Y$}\right.}\right),& {{\displaystyle \sum }}_{n=1}^{N}Pn \left(i,j\right)>Y\hfill \\ \hfill ln\left({{\displaystyle \sum }}_{n=1}^N{e}^{\raisebox{1ex}{$P_n\left(i,j\right)$}\!\left/ \!\raisebox{-1ex}{$Y$}\right.}\right),& {{\displaystyle \sum }}_{n=1}^{N}Pn \left(i,j\right)<Y\hfill \end{array}\right.$$

(3)

where the parameter $\alpha$ is the intensity coefficient of the enhanced image. The value of $\alpha$ is set to 1.35. The parameter $Y$ is the threshold to the Hu of the soft tissue. $Pn\left(i,j\right)$ is the $i$ row $j$column value of the $n$th image in the curved MPR image set. $N$ is the number of images.

$$P=\left\{\begin{array}{cc}\beta \ast G(P_1),& P_0<S\hfill \\ \gamma \ast G\left(P_0\right)+\left(1-\gamma \right)\ast W(P_0),& P_0>S\hfill \end{array}\right.$$

(4)

where the parameter $\beta$ is an adjustment factor that controls the contrast of the non-area of interest, $\beta =1.3$. The parameter $S$ is the area of interest threshold in the enhanced image, $S=105$. The parameter$\gamma$ is the weight factor used to control the enhanced details, $\gamma =0.4$.$W$ represents the bilateral filter function. $G$represents the Gaussian filter function.

After the generation of$P_0$, the region of interest in the enhanced image is mapped to the panoramic image in the new intensity interval by Formula (4) to improve the contrast of different regions of the image. The pseudo-code of the image enhancement and synthesis algorithm is shown in Algorithm 1. Lines 1–8 and lines 9–14 in Figure 6 represent Formula (4) and Formula (5), respectively.

Algorithm 1 Image synthesis and enhancement algorithms

Input: MPRSets(Psets)

Output:P

1: N←Num(Psets)

2: for Pn(i, j) ∈ Pn do

3:     if Pn=1→N(i,j) > Y then

4:        P0(i,j) = αYln(ePn=1→N(i,j)/Y)

5:     else

6:       P0(i,j) = ln(ePn=1→N(i,j)/Y)

7:     end if

8: end for

9: if P0(i,j) > Y then

10:   P(i,j) = βG(P1)

11: else

12:    P(i,j) = γG(P0) + (1-γ)W(P0)

13: end if

14: return P

3. Results and Discussion

3.1. Experimental Environment and Implementation Time

The experiments in this paper were all performed on a computer with Intel(R) Xeon(R) E-2136 [email protected] GH (Dell, Round Rock, TX, USA) and 32GB RAM (Kingston, Fountain Valley, CA, USA), using Python 3.8.5 (G. van Rossum., Amsterdam, The Netherlands). We performed the image reconstruction experiment and recorded the time after the CBCT data reading was completed. The average time spent on the method based on the reconstruction of the dental arch area was 4.58 s. The method proposed in this paper has a maximum time of 5 s and an average calculation time of 3.04 s. This shows that the reconstruction method based on multi-arch curves can effectively reduce the computational complexity and improve efficiency.

3.2. Dental Arch Detection Effect

There are several 3D imaging software currently available for reconstructing panoramic images, such as SimPlant (Materialise Inc. Rotherham, UK) and eXam Vision (KaVo Dental GmbH, Biberach, Germany). In eXam Vision, the software supports the automatic creation of dental arch curves or manual creation of dental arch curves, and a panoramic image corresponding to the curve can be obtained. Therefore, we used the same CBCT data to compare the panoramic images corresponding to the dental arch curves generated by the eXam Vision software and our proposed dental arch detection method. It can be seen in Figure 6a–c that the panoramic image reconstructed by the software failed to show some incisors and molars. Figure 6d–f are panoramic images generated by the dental arch detection method in this paper, and the pictures show the entire dentition more completely. We believe that the key factor for improving the accuracy of the proposed dental arch detection method over different ranges regarding the performance of dental arch curve fitting is to apply sub-regional detection for reducing disturbance among various ranges and add a control point adjustment panel for decreasing fitting error.

3.3. Different Types of Panoramic Image Effects

In addition, when the patient has implants such as dental crowns in the mouth, metal implants can compress the brightness distribution of other tissues on the image, resulting in an overall darkening of the final panoramic image. This is due to the strong ability of the metal to absorb X-rays and the weak absorption ability among tissues, so we increased the strength of the tissue areas outside of the metal species plants, redistributing them to new areas, and we took the enhanced image that improves the clarity of the boundary as the foreground and the original panoramic image after Gaussian blurring as the background to increase the depth and spatial contrast of the final panoramic image.

The images of experimental results for common oral CBCT data are presented in Figure 7. These CBCT data included normal dentition, missing teeth, restorations with crowns, implants, and mixed cases, which covered most types of dental implant cases in different locations. In these panoramic images, the tooth area and the treatment area have better luminance performance. Figure 7a show a panoramic image of normal oral data generated with good contrast. Panoramic images generated from CBCT data of implant and filler treatments are shown in Figure 7b–f,i, the orientation and depth information of the maxillary and maxillary implants can be observed, and the treatment area can be clearly distinguished from the surrounding tissue with clear shapes, and distinct layers can be seen in the image. Panoramic images of crown restoration are shown in Figure 7e–i. The panoramic image does not appear black due to the high brightness of the restored crown, maintaining a balanced contrast between the treatment area, the tooth area, and the surrounding tissue, and does not require the user to adjust the window width and window level to observe the oral tissue information. It can be seen that our proposed method has better brightness and contrast performance in the tooth area and the treatment area in various treatment situations.

3.4. Quality Evaluation of the Panoramic Radiographs

In this experiment, we divided the image into three regions: irrelevant tissue region, tooth region, and treatment region. By referring to the quality requirements of CBCT images [27,33,34,35], the evaluation content of panoramic images is determined (

Table 1. Panoramic image score descriptions.

Score	Evaluation Content
1	The differences in brightness and contrast between different regions are too large, and the structures within each region are not visible.
2	The contrast difference between different regions is not obvious, and the boundaries of each region are not obvious.
3	Brightness and contrast between different tissues are good.
4	There is good density and contrast between different tissues. Brightness is uniform across areas, and internal details are visible.

). The evaluation results were scored by three dentists with extensive clinical experience based on the evaluation content, with a score ranging from 1 to 4 points, with a full score of 12 points. An average score was calculated to represent the contrast of each panorama image.

Panorama images generated from four different methods were evaluated by three dentists. The rating process is conducted using blinded method. The results show that the panorama image reconstructed by our proposed method has the highest overall score of 11.16 ± 2.64. The scores assessed by the dentists are shown in Figure 8.

3.5. Limitations

Although our proposed method achieves excellent performance on the task of automatically reconstructing panoramic images, it also has some limitations. First, only a small amount of curved MPR images in the dental arch region were extracted for the growth shape of the tooth; thus, the internal continuous information collected was insufficient. Additional information on the growth shape of the tooth could potentially benefit the display and further improve the quality of panoramic images. Second, the values of α, β, and γ in Formulas (3) and (4) may vary due to different models of CBCT scanners, consequently increasing the workload. Finally, we have only tested the algorithm on limited data and have not yet applied it clinically.

4. Conclusions

In this paper, we proposed a fast, automatic reconstruction method for panoramic images based on CBCT. This method first detects the dental arch curve through axial MIP images over various ranges and then applies an enhanced image fusion algorithm to extract different tissues and adjust the intensity distribution during the reconstruction process. Our results indicate that the proposed method effectively reduces the disturbance among different types of tissue by utilizing sub-regional detection, notably increases the fitting accuracy of the dental arch curve by implementing a control point adjustment panel, and significantly brings down the influence of metal implants through intensity distribution adjustment. Therefore, the dentition display appears on ideal brightness, high level of completion, and proper contrast in both tooth and treatment areas. Despite the limitations mentioned in Section 3.5, our work presented a high-speed, reliable, and satisfactory panoramic image reconstruction technique for digital dentistry.

In future work, we plan to construct new transfer functions for different tissues or regions of interest during the synthesis of panoramic images or to use different algorithms to enhance and improve panoramic images [36,37,38]. In addition, we will try to extract the dental arch curves at different positions to further improve the performance and generality of the algorithm.