A Novel Thresholding Based Algorithm for Detection of Vertical Root Fracture in Nonendodontically Treated Premolar Teeth.

In this paper, an efficient algorithm is proposed for detection of vertical root fractures (VRFs) in periapical (PA), and cone-beam computed tomography (CBCT) radiographs of nonendodontically treated premolar teeth. PA and CBCT images are divided into some sub-categories based on the fracture space between the two fragments as small, medium, and large for PAs and large for CBCTs. These graphics are first denoised using the combination of block matching 3-D filtering, and principle component analysis model. Then, we proposed an adaptive thresholding algorithm based on the modified Wellner model to segment the fracture and canal. Finally, VRFs are identified with a high accuracy through applying continuous wavelet transform on the segmented radiographs and choosing the most optimal value for sub-images based on the lowest interclass variance. Performance of the proposed algorithm is evaluated utilizing the different tested criteria. Results illustrate that the range of specificity deviations for PA and CBCT radiographs are 99.69 ± 0.22 and 99.02 ± 0.77, respectively. Furthermore, the sensitivity changes from 61.90 to 77.39 in the case of PA and from 79.54 to 100 in the case of CBCT. Based on our statistical evaluation, the CBCT imaging has the better performance in comparison with PA ones, so this technique could be a useful tool for clinical applications in determining the VRFs.

INTRODUCTION

Vertical root fractures (VRFs) are kind of fractures, which transverse the length of tooth from the crown to the apex, usually affect one side of the root, and appear on the facial and lingual surfaces of anterior and posterior teeth.[1] VRFs occur in all teeth due to a centrifugal force caused by factors such as occlusion, trauma, main restorations, the additional pressure during endodontic therapy, and weak positions.[2] Most often VRFs are encountered in endodontically treated teeth; however, it may occur in intact teeth. The prevalence of VRFs in endodontically treated teeth is much higher than in nonendodontically ones (vary from 2% to 5%).[3] Early diagnosis of root fracture can prevent from the spreading the damage to the surrounding tissues such as probing defect, apical bone loss, localized widening of periodontal ligament (PDL), pulpal lesion, and inflammation.[4]

In general, there are two important parameters in the identification of VRFs; clinical signs and radiographs. VRFs in periapical (PA) radiography are diagnosed by the presence of isolated pockets at one side of the tooth, looseness in the fibers of the PDL, additional space in the vicinity of root canal filling, cement protrusion in the canal, and bone loss at one side of the tooth. However, in this type of radiographs, fracture can be observed only when X-rays pass the fracture line. Thus, we generally need two or three radiographs with different angles. According to a study conducted by Nair et al.,[4] the fracture can be detected when X-rays meet the fracture line at an angle of 4°.

Several studies have reported successful results on the detection of VRFs by utilizing high-resolution imaging systems such as tuned aperture computed tomography (CT), optical coherence tomography, and local CT.[5] In spite of their advantages, some of these systems were not commercially available,[4] or they were designed only for medical purposes that were not suitable for dental use. Furthermore, unsatisfactory image quality with lower resolution was reported when compared with cone beam CT (CBCT).[67]

Image processing algorithms could be utilized as a powerful tool for the diagnosis of VRFs in PA and CBCT radiographs. Kositbowornchai et al.[8] proposed a probabilistic neural network model for diagnosis of VRFs in nonroot treated teeth. In this method, four horizontal lines crossing the root in the radiographs of healthy and fractured teeth were considered. One of the key challenges of this algorithm is that, not only it examines the fractures for untreated teeth, but also it has limitations on finding the location of them. Hence, this study was investigated only the premolar teeth with a single root and was not focused on endodontic filling materials or radiopaque post.

By highlighting the importance of the diagnosis of VRFs and the limitations of PA radiographs for detection of VRFs, this paper proposes an efficient algorithm based on continuous wavelet transform (CWT) and adaptive thresholding for detection of VRFs in untreated teeth of CBCT and PA radiographs. The rest of the paper is organized as follows; Section 2 elaborates on procedures of teeth preparation and the imaging protocol. The proposed algorithm and its various steps are described in details in Section 3. Results of the algorithm implementation in the designed user-friendly software environment and the discussion of the proposed algorithm based on evaluation criteria are presented in Section 4. Finally, the conclusion of the paper is provided in Section 5.

TEETH PREPARATION PROCEDURES AND IMAGING PROTOCOL

A total of 40 radiographs of noncarious nonroot treated premolar teeth, which were kept in formalin for consolidation after cleaning, were examined by a microscope with a ×20 magnification to ensure the absence of fracture. Teeth crown was cut-off from 2 mm below the cementoenamel junction. Canal preparation was conducted using the rotary system named ProTaper up to number F5, and the filling was carried out using F5 ProTaper Gutta-percha points. The teeth were fixed in an acrylic resin box using 1 mm of wax. Then, teeth were attached to the universal testing machine with a 60° bevel and were exposed under pressure to achieve fracture. After achieving fracture, the teeth were taken out and re-examined by a microscope with a ×20 magnification, this time, to ensure the presence of fracture.

RVG sensors were used for PA radiography. Images were taken with the teeth being faciolingually in parallel to the sensor and were saved in JPG format – resolution: 1320 (horizontal/vertical) × 1024 (depth) × 1536 (height). X-ray tube was emitted in the form of a round collimator with a maximum voltage of 65 Kv and a maximum current of 8 mA for 0.08 s. The distance between the teeth and the tube was 24 cm and the distance between the object, and the film was 1 cm. In CBCT imaging, slices were made at axial levels with a distance of 1 mm to achieve threshold equalization. The CBCT system used a cone-shaped X-ray beam and was equipped with a flat-panel detector (pixel array 1536 × 1920), 360° rotation, 18s scanner, and a maximum voltage of 110 Kv. The initial and final reconstruction was carried out by the NNT Viewer Software, version 2.17 (http://www.newtom.it/en/products/newtom-nnt), and exposure conditions of the device were set automatically.

PA/CBCT radiographies were divided into three categories as small, medium, and large based on the fracture space between the two fragments. Table 1 shows the different situations of this division.

Table 1

Classification of periapical/cone beam computed tomography graphies into three sub-categories based on the fracture space between the two fragments

METHOD

Figure 1 shows a block diagram of the proposed algorithm. The main steps of the algorithm are discussed in detail as follows.

Figure 1

Block diagram of the proposed algorithm

Image Denoising Using a Hybrid Block Matching 3-D Filtering - Principle Component Analysis Model

The proposed algorithm uses the block matching and 3-D filtering and principle component analysis (PCA) hybrid model for denoising the image.[9] The model is summarized in Figure 2. Let assume that the input image is destructed by additive white Gaussian noise with zero average and σ2 variance.[9] The input image is scanned as an array. Then, the following procedures are applied to each processed pixel:

Figure 2

Block diagram of denoising in the proposed algorithm

For each processed pixel, a shape adaptive neighborhood can be found with the original pixel as the center of the neighborhood using the 8-connected LPA-ICI model based on Lee et al.[9] The neighborhood is located within a square block of a fixed size. This block is called the reference block. The number of pixels in the neighborhood is represented by Nel

Any block similar to the reference block can be found using block-matching. The shape adaptive neighborhood is extracted from among these matched blocks using the method described in the previous step. The number of matched blocks is represented by Ngr

Determining a transform and applying it to the neighborhoods of an adaptive shape. To this end, we consider a thresholding level (τ) and investigate two conditions:

First condition: If we have , it means that an acceptable number of matched neighboring pixels are selected for estimating the PCA matrix. This matrix forms the basic vectors of the PCA model. So, only eigenvectors whose corresponding eigenvalues are greater than the defined threshold are selected at this stage

Second condition: If we have , it means that an acceptable number of matched neighboring pixels are not selected for estimating the PCA model. So, a specific eigenvalue is selected for the eigenvector as proposed by Lee et al.[9]

Forming a three-dimensional array by connecting shape adaptive neighborhoods (min [Ngr, N2]) to a reference block which is most similar. N2 is a constant parameter that limits the number of filtered neighborhoods

Applying transform in step 3 to each group of shape adaptive neighborhoods. At this stage, a one-dimensional orthogonal transform (i.e., wavelet transform) is applied to each three-dimensional group

Applying a hard-thresholding to the three-dimensional group for achieving image shrinkage

Applying three-dimensional inverse transform from step 5 to find the estimates in each group of shape adaptive neighborhoods, and

Relocation of collected estimates back into their original positions using weighted averaging.

Image Segmentation for Determining the Fractures Using the Adaptive Thresholding Algorithm

VRFs are considered as the edges of the image. So in this paper, we propose an algorithm for segmentation based on adaptive thresholding for determining the fractures in CBCT and PA radiographs. The idea of segmentation is an extension of the Wellner method[10] in which any pixel is compared to its adjacent pixels. Thus, if the value of the denoised pixel is T percent below the Wellner method's output value, it belongs to the black pixel. Otherwise, it belongs to the white pixel. The main problem of the Wellner method in PA and CBCT radiographs is its dependence on the sweeping up of pixels. In addition, in this type of image, this method alone cannot be an appropriate criterion for displaying the existing pixels in the image neighborhood, since all neighborhood samples are not dispersed in all directions with equal distribution.

To overcome the Wellner algorithm's issue with PA and CBCT radiographs, we have conducted two important modifications. At first, a maximum filter is applied to the denoised image. The output of the maximum filter is the maximum value of a pixel in its neighborhood in a window. The dimension of the selected window is set by 10 × 10. The second modification is applied to eliminate the dependence of the thresholding method on the sweeping up of pixels, in which the whole image is integrated using Eq. 1:

where, I (x, y) is the input image and f (x, y) is the result of input image integration. Therefore, instead of calculating the average of the last s pixel, the maximum value is first calculated in a window with each pixel as the center. Then a window with w × w dimensions is set out on each pixel obtained from the image for calculating its average value.

The maximum filter highlights the effects of the edge in the image. Edges include canal, fractures, and the complete configuration of the teeth. In PA radiographs, canal is produced by a lower level of intensity transition compared to fracture. This feature causes the fractures to be more clearly visible than internal teeth borders (i.e., canal) by changing the threshold level in segmentation. Figures 3-6 shows the results of segmentation using the proposed algorithm and its comparison with the conventional Wellner algorithm for two samples of PA radiographs (with medium and small fractures) at different thresholding values (sample 1: T = 0.89, T = 0.79 and T = 0.69; sample 2: T = 0.95, T = 0.80 and T = 0.71). By comparing Figures 3 and 4 as well as 5 and 6, we can see the effect of the maximum filter in the Wellner method. In the proposed algorithm, the fracture is visible in Figure 4c (for T = 0.69), but the conventional Wellner algorithm fails to identify the fracture under similar conditions. This becomes even more important in the diagnosis of small fractures.

Figure 3

Results of segmentation of periapical radiographs with medium fractures via the conventional Wellner algorithm at different thresholding values. (a) T = 0.89 (b) T = 0.79 and (c) T = 0.69

Figure 6

Results of segmentation of periapical radiographs via the proposed algorithm at different thresholding values. (a) T = 0.95 (b) T = 0.80 and (c) T = 0.71

Figure 4

Results of segmentation of periapical radiographs via the proposed algorithm at different thresholding values. (a) T = 0.89 (b) T = 0.79 and (c) T = 0.69

Figure 5

Results of segmentation of periapical radiographs with small fractures via the conventional Wellner algorithm at different thresholding values. (a) T = 0.95 (b) T = 0.80 and (c) T = 0.71

The variance of intensity that forms edges in CBCT images is the same in both canal and fractures and the maximum filter highlights the edges as well. Results of segmentation for a CBCT image at different thresholding values are presented in Figures 7 and 8. As can be seen, the edge to surface intensity deviations is high in this type of imaging, and this leads to more highlight the edges.

Figure 7

Results of segmentation of cone beam computed tomography images via the Wellner algorithm at different thresholding values (a) T = 0.87 (b) T = 0.77 and (c) T = 0.57

Figure 8

Results of segmentation of cone beam computed tomography images via the proposed algorithm at different thresholding values (a) T = 0.87 (b) T = 0.77 and (c) T = 0.57

Adaptive Thresholding Algorithm

To achieve the best segmentation, we need the best threshold image. The best threshold image is a coefficient of the maximum/average filter's output in PA/CBCT radiographs. To determine an efficient coefficient to detect the threshold image, we propose an adaptive thresholding algorithm. This algorithm is based on a comparison between the reference image and the image obtained from different thresholding values. Steps of the proposed adaptive thresholding algorithm are summarized as follows.

Start

Acquiring threshold and segmented images for a fixed T

Applying CWT to images from step 1 and creating image masks

Comparing the output image with the reference image and calculating the interclass variance of the two images based on the Eq. 2:[11]

where, μT is the average of the entire image; μ1(t) and μ2(t) are the averages of the first and second classes, respectively; and w1(t) and w2(t) are the weights of each class, respectively.

Repeating the above steps for different T values, until the minimum interclass variance is obtained.

End

Output: Efficient T value and the number of sub-images.

Feature Extraction by Applying Continuous Wavelet Transform

Wavelet transform is used as a powerful tool in various signal and image processing functions. Two-dimensional wavelet transform breaks a signal down to a set of basic vectors all of which are scaled, shifted, and rotated versions of mother wavelet functions.[1213] Assume g() and, = (x,y) representing pixel intensity and image coordinates, respectively. CWT is defined as follows:[14]

where = (b,b) represents the transition parameter showing wavelet shift in the spatial domain; Φ represents the rotation parameter.[13] With the rotation function defined as r−Φ(x, y), wavelets rotate in the angular domain at a Φ degree angle. So we have:

in Eq. 3 a is the scale parameter which corresponds with the spatial frequency of wavelet functions. ψ* represents the complex conjugate of mother wavelet (ψ) functions. (Cψ) represents a normalized constant such that the Fourier transform of the ψ function meets the acceptability condition as Eq. 5:

Where FT (ψ(·))represents the Fourier transform of the ψ function and = (,k) represents the spatial frequency in a nondimensional R2 space. The proposed algorithm used the esmexh mother wavelet for calculating wavelet transform coefficients. The esmexh mother wavelet is defined in the frequency and spatial domain as follows:[14]

Where σ is the variance of the wavelet from its average value and T is expressed as follows:

Where ε is a constant value and is set by 0.5. The reason for using the esmexh mother wavelet is that the features extracted from the image, which includes the edges of the image, are part of singular points. Therefore, the mother wavelet that has a greater number of zero moments is more conducive to the highlighting of the edges. In the mother wavelet function, four scales and 16 different phases were used for PA radiographs. Also for reconstruction, between 5 and 10 sub-images were gathered together and the efficient value of sub-images was determined based on the adaptive thresholding method. In CBCT, four scales are used. However, the sub-images were considered in eight phases and between four and eight sub-images were gathered together. The reason of choosing the necessary angles to obtain the sub-images lies in the resolution of fractures in any radiograph. Hence, in PA images where it is more difficult to diagnose fractures and the internal edges (fractures) are thinner, we need smaller angles for highlighting. However, the exact opposite applies to CBCT radiographs. Figure 9 shows the results of CWT on CBCT and PA radiographs.

Figure 9

Results of continuous wavelet transform on periapical and cone beam computed tomography radiographs

Locating the Vertical Root Fractures

After applying the CWT, the reconstructed image is masked and transformed into a binary structure for identifying the exact location of the fracture. The basis for the binary transformation of the image is thresholding, such that all pixels whose intensity exceeds the threshold value (i.e., 0.21 in this the paper) will be equal to one and the rest will be equal to zero. To produce the mask, the resulting binary image is multiplied by the image of the wavelet. Figure 10a and b show the results of applying the mask and the final result in CBCT and PA radiographs, respectively.

Figure 10

The mask produced for identifying the exact location of the fracture and the final result in (a) periapical radiographs, and (b) cone beam computed tomography radiographs

RESULTS AND DISCUSSION

Figures 11 and 12 show the software package designed for detecting the location of fractures in CBCT and PA radiographs. This software package was designed and implemented in the MATLAB software (http://www.mathworks.com/products/matlab/) environment by the Department of Oral and Maxillofacial Radiology, the Faculty of Dentistry, Tabriz University of Medical Sciences. Furthermore, reference images were produced by radiologists and were used for comparison with the results of the proposed algorithm.

Figure 11

The software package designed for detecting the location of fractures in periapical radiographs

Figure 12

The software package designed for detecting the location of fractures in cone beam computed tomography radiographs

According to the classification stated in Table 1, in PA radiographs, the fractures in teeth are classified into small, medium, and large. Besides, all cases of CBCT images in this paper have the fracture space higher than 2.5 pixels; hence, these radiographs are classified as large fracture. Based on this explanation, in the graphies shown in Figure 13, PA1, PA2, PA3, and PA4 are small, the fractures in teeth PA5 and PA6 are medium, and the fractures in teeth PA7 and PA8 are large. The final results of our proposed algorithm are presented in Figure 14. As can be seen, the proposed algorithm can detect the exact location of small, medium, and large fractures in all radiographs. It should be noted that in radiograph PA1, the horizontal line relates to image artifacts and is not a fracture. Besides, radiograph PA8 shows a superimposed fracture in the buccal or lingual surface and the tooth is not filled with Gutta-percha.

Figure 13

Samples of periapical radiographs with vertical root fracture (PA1, PA2, PA3 and PA4 with small vertical root fracture; PA5 and PA6 with medium vertical root fracture, PA7 and PA8 with large vertical root fracture)

Figure 14

Final results of applying the proposed algorithm on periapical radiographs with vertical root fracture (PA1, PA2, PA3 and PA4 with small vertical root fracture; PA5 and PA6 with medium vertical root fracture, PA7 and PA8 with large vertical root fracture)

In Figures 15 and 16, CBCT radiographs and the final results of the proposed algorithm are presented, respectively. These radiographs are labeled by CT1 to CT7.

Figure 15

Samples of axial cone beam computed tomography radiographs with vertical root fracture

Figure 16

Final results of applying the proposed algorithm on axial cone beam computed tomography radiographs with vertical root fracture

To evaluate the performance of the proposed algorithm we have used the evaluation criteria as follows:

Sensitivity (Sn), specificity (Sp), precision (P), and approximation correlation (AC).

These parameters are defined as follows:[15]

where TP and TN represent the number of pixels correctly identified as fracture and nonfracture regions, respectively. Similarly, FP and FN represent the number of pixels mistakenly identified as fracture and nonfracture regions, respectively. We have utilized the above evaluation criteria to assess the correct detection rate.

According to Eqs. 8-11, the Sn parameter provides a measure of fracture pixels correctly identified as fracture regions. The Sp parameter also provides a measure of fracture pixels correctly detected by the proposed algorithm as fracture regions. Accordingly, both these parameter (Sn and Sp) can be considered in terms of conditional probability. Finally, the precision parameter shows the accuracy of the segmentation algorithm.

Negative predictive value (NPV), F_measure and G_mean parameters are defined in Eqs. 12-14:[16]

By considering P and NPV, F_measure examines the precision of the algorithm and assigns it a score– one indicates the highest and zero represents the lowest score. This parameter is used in many cases including border detection, segmentation, and clustering.

Peak signal to noise ratio (PSNR): Since many signals have a wide dynamic range, this measure is expressed logarithmically. This measure is often used for measuring the quality of images after reconstruction. A higher PSNR thus means a better-reconstructed image quality. If the noiseless image is represented by I and the noised image by K, we have:[1718]

where MAXI represents the maximum possible pixel value in the image. If there are eight bits per pixel, the maximum pixel value is 255. Otherwise, it is calculated in the following manner:

MAX = 2 –1      (16)

Receiver operating characteristic (ROC) Curve: ROC curve was developed in 1950 as a technique for organizing classifiers and visualizing their performance.[19] In determining fracture zones, the ROC curve evaluates the impact of TP and FP by considering an arbitrary thresholding level and is defined as an XY-diagram of the TP rate as a function of the FP rate at different thresholding level values. The area under the ROC curve is a measure of the performance of fracture/nonfracture separation system. A greater area under the ROC curve in a certain algorithm represents a higher accuracy. It should be noted that trapezoidal integration was used for calculating the area under the ROC curve. This measure is equal to the fact that a classification is more likely to rate a positive random value than a negative random value.[20]

As can be seen in Table 2, a σ = 15 yields the best result for image denoising in PA radiographs such that in the case of the PA4 radiograph, which is considered as an short VRF, the PSNR value is improved for 9.2 and 8.36 at σ = 15 in comparison with σ = 5 and σ = 25, respectively. This superiority can also be seen in the case of other radiographs as shown in Table 2. It should be noted that in the PA5 radiograph, which is considered as medium VRF fracture, the PSNR value is obtained as greater at σ = 25 in comparison with σ = 15. However, the latter sigma value (σ = 15) is selected due to an equal Sp value at these two sigma values (Sp = 99.92). In CBCT radiographs, which are considered as large VRF, as presented in Table 3, the best sigma value is different for each image. However, with a good approximation, it can be concluded that σ = 15 is the best match. In this study, the aim of denoising is to sharpen the surfaces and consequently improve the accuracy of edge detection. However, if denoising is carried out with low sigma values, the surfaces will not become sharp. High-sigma values also destroy the edges and reduce diagnostic precision. By comparing Tables 2 and 3, we find that changes in the sigma vale have a greater impact on the PSNR value in PA radiographs due to the presence of medium and small fractures. In contrast, CBCT radiographs are characterized by a greater stability in the face of noise when compared with PA radiographs. It should be noted that the results presented in Tables 2 and 3 have been achieved based on the best thresholding level, which was obtained by changing the image threshold coefficient (T) (from 3.0 to 9.0).

Table 2

Quantitative values of peak signal to noise ratio, precision and specificity in periapical radiographs at various sigma values and the best image threshold coefficient

Table 3

Quantitative values of peak signal to noise ratio, precision and specificity in cone beam computed tomography radiographs at various sigma values and the best image threshold coefficient

The quantitative values of evaluation parameters for CBCT and PA radiographs are presented in Tables 4 and 5 for the best sigma value. In addition, ROC curves for PA and CBCT radiographs are presented in Figures 17 and 18, respectively. In PA radiographs, the area under the ROC curve for large fractures is improved for 1.02 and 1.1 in comparison with medium and small fractures, respectively. The superiority of the proposed algorithm is also visible in CBCT radiographs such that the area under the ROC curve is equal to 0.9755 in the CT5 radiograph and equal to 0.9607 and 0.9385 in CT1 and CT6 radiographs, respectively. From Figures 17 and 18, it can be concluded that the proposed algorithm has a better performance in CBCT compared with PA radiographs. This is also visible in Tables 4 and 5 such that the Sn parameter in CBCT radiographs is higher compared with PA radiographs for the best image threshold coefficient. In PA radiographs, the Sn parameter ranges from 61.90 (for PA2) to 77.39 (for PA6). However, the Sn range is between 79.54 (for CT6) and 100 (for CT3 and CT7) in CBCT radiographs. Figures 19 and 20 present the Sp parameter as a function of the thresholding level in PA and CBCT radiographs. As can be seen, Sp values are almost constant for all thresholding levels in the proposed algorithm. The range of Sp is between 99.47 (for PA8) and 99.92 (for PA3 and PA5) in PA radiographs and between 98.25 (for CT5) and 99.81 (CT3). This shows that the proposed algorithm has a good stability for locating VRFs in clinical applications.

Table 4

Quantitative values of evaluation parameters in periapical radiographs for the best thresholding level

Table 5

Quantitative values of evaluation parameters in cone beam computed tomography radiographs for the best thresholding level

Figure 17

Receiver operating characteristic curve for periapical radiographs

Figure 18

Receiver operating characteristic curve for cone beam computed tomography radiographs

Figure 19

Sp curve as a function of thresholding level for periapical radiographs

Figure 20

Sp curve as a function of thresholding level for cone beam computed tomography radiographs

CONCLUSION

A novel thresholding based algorithm for segmentation of VRFs in PA and CBCT radiographs in nonendodontically treated premolar teeth has been proposed. Results show that the proposed algorithm has good stability in finding the location of VRFs and can accurately locate any types of fractures – either small, medium, or large – in CBCT and PA radiographs. As it can be found from the results, the range of Sp deviations for PA and CBCT radiographs are from 99.47–99.92 to 98.25–99.8, respectively. This proves that our algorithm could be a useful tool for clinical applications and also could be utilized to accurately localize the VRFs. One of our limitations in this work is the use of premolar teeth, with a single root and no endodontic filling material or radiopaque post. Our goal in future researches is to study multi-rooted teeth for detection of VRFs.

Financial Support and Sponsorship

Nil.

Conflicts of Interest

There are no conflicts of interest.