Measure the Dental Chamber Volume with DICOM Images from Cone-Beam Computed Tomography Can Be Improved with a Simple Algorithm

1. Introduction

Cone beam computed tomography (CB-CT) has been introduced succesfully in daily clinical practice. Small-field CB-CT allows to obtain teeth three-dimensional images of good quality without overlapping neighboring anatomical structures. This, associated with low radiation levels, has made CB-CT an image technique routinely used in diagnosis and planning of dental treatments ([1,2,3]). One of the applications of CBCT is the measurement of the pulp chamber from DICOM images. However, the values acquired show slight discrepancies with those obtained by micro-computed tomography (micro-CT). Micro-CT is currently regarded as the gold standard technique for assessing pulp chamber volume ([4,5,6,7,8]). In addition, several factors limit the use of micro-CT in daily dental practice, including the high price of data analysis equipment and software, the time required for digital scanning and reconstruction, the volume limitations it can analyze, and the high doses of radiation needed ([1,5,9]).

Therefore, the availability of reliable tools to measure the volume of the pulp chamber from CBCT images, obtaining values as close as possible to those obtained by micro-CT, could be of great interest in daily clinical practice. Some studies are available that corroborate its accuracy and reproducibility de estos estudios comprativos ([10,11]).

Therefore, the aim of this work was to design a simple algorithm capable of calculating the volume of the pulp chamber from DICOM images obtained with CB-CT. The results obtained show that its application allows to obtain values almost identical to the real values of micro-CT. This is a validation of a method, not a quantitative study of a large series of samples.

This study was designed to validate volume of the pulp chamber obtained using CB-CT in vivo (PromaxR 3D Max CBCT) by comparing it with those obtained using micro-CT ex vivo. Dental treatments involving dental chamber are becoming personalized and the same tooth differs from one patient to another ([12]). Thus, the exact knowledge of the anatomy a volume of the dental chamber is essential to achieve successful treatments.

2. Materials and Methods

Five permanent molars, corresponding to 5 subjects (3 females and 2 males), with ages ranging between 59 and 66 years, were analyzed. Teeth used in the study were programmed for exodontia due to different odontological treatments and were free of pulp camber damage. This is the reason why the pulp was intact and free of damage despite the age of the patients. It should be considered that teeth with caries suffer retraction of the pulp chamber, decreasing the volume of the pulp chamber. Likewise, restorations of decayed teeth produce artifacts in the images that prevent the volume from being correctly assessed.

Pieces were collected under written informed consent and obtained in accordance with Spanish legislation (RD 1301/2006; Law 14/2007; RD 1716/2011; Order ECC/1404/2013). Before exodontia CB-CT images were obtained using a tomographic equipment ProMax3DMid ProFace, with a voltage of 120kV and a tube current of 5mA, except for the tooth identified as 5, which were 90kV and 10mA. For the analysis of DICOM (digital imaging and communication in medicine) images, the software used was Planmeca Romexis 6.2.1.25. Briefly, measurements were made as it follows: the axial and transverse axes were oriented according to the direction of the tooth. The tool "Measure a cube" was used to trace the area that covered the pulp chamber in the 3 visions: axial, panoramic and sagittal/orthoradial. Selecting a point inside the cube belonging to the pulp chamber as the reference density for the "3D Growth Region" tool with the preset "Root Cavity" adjusting the threshold favorably (Figure 1 a-d).

After exodontia teeth were preserved in physiological saline solution at 4o C until use for micro-CT study, which was carried out at Centro Nacional de Investigacion sobre la Evolucion Humana (CENIEH; Burgos, Spain). The images were obtained with V | Tome | X s 240 de GE Sensing Inspection Technologies Phoenix X-Ray and for the analysis of the images, the software Dragonfly was used. The pulp canal was deleted from the rest of the tooth by adjusting the histogram. Once the pulp canals were isolated, a box was defined that corresponded to the region of the pulp chamber (Figure 1e-i). The software itself has a function that directly provides the volume value of the selected region (connected components; Figure 1i).

To calculate the pulp chamber volume, 300 DICOM images for each of the 5 samples, corresponding to different axial planes, were available. Nevertheless, only the images corresponding to the pulp chamber were selected, resulting in a reduction of 14/300 for sample 1, 16/300 for sample 2, 9/300 for sample 3, 16/300 for sample 4, and 12/300 for sample 5. The selection was based on a meticulous visual analysis of each image by two experts, who established the initial and final planes containing parts of the pulp chamber. The image chosen as the initial one was the closest to the pulp canals among those that contain the pulp chamber, and the rest were selected in an ascending order. On the other hand, the voxel size varied for each sample, with a side length of 0.25 mm for the first three samples, 0.2 mm for the fourth sample, and 0.15 mm for the fifth sample.

3. Design of the algorithm

The effectiveness of image segmentation algorithms can vary depending on a range of factors, such as image complexity, image quality, the type of features being sought, noise present in the image, and the required accuracy for the specific task. In this case, a simple algorithm has been implemented to detect pixels belonging to the pulp chamber. Like any other image segmentation algorithm, it is based on the discontinuity or similarity between the gray levels of neighboring pixels, allowing it to locate changes in the gray level of adjacent pixels. This algorithm can find isolated pixels, detect edges, and divide the image into zones that have similar gray values based on a certain threshold.

Figure 2. CB-CT images corresponding to the molar identified as sample 2. Axial images and area occupied by the pulp chamber (highlighted with red dots). In this sample the algorithm values match those of micro-CT.

Figure 2. CB-CT images corresponding to the molar identified as sample 2. Axial images and area occupied by the pulp chamber (highlighted with red dots). In this sample the algorithm values match those of micro-CT.

3.1. Preliminaries

The images used in this study are matrices in which each element corresponds to a position in the image (pixel). The value of each matrix element is determined by the color of the corresponding pixel. Since the images being used are in gray scale, each image is determined by a single matrix. Once the images are uploaded, an automatic check is performed to determine if the image has borders, and if so, they are not considered in the execution of the algorithm. This way, unnecessary calculations are avoided without preprocessing the images in the dataset.

Figure 3. CB-CT images corresponding to the molar identified as sample 5. Axial images and area occupied by the pulp chamber (highlighted with red dots). In this sample the algorithm values were very different from those of micro-CT.

Figure 3. CB-CT images corresponding to the molar identified as sample 5. Axial images and area occupied by the pulp chamber (highlighted with red dots). In this sample the algorithm values were very different from those of micro-CT.

Defining the application d as follows:

$$\begin{gathered} \mathrm{d} : {\mathit{Z}}^2\to \mathit{R} \\ \mathit{x},\mathit{y} \to d(\mathit{x},\mathit{y})=\left|x_1-y_1\right|+\left|x_2-y_2\right| \end{gathered}$$

(1)

where $\mathit{x}$ = ($x_1,x_2$), $\mathit{y}$ = ($y_1,y$) ∈ ${\mathit{Z}}^2$, then (${\mathit{Z}}^2$, $\mathrm{d}$) is a metric space.

Given α ∈ N, on a set A ∈ $M_{2\times N}$(Z), the equivalence relation ${~}_{\alpha }$ is defined:

$$\mathit{x},\mathit{z}\in A \underset{\square }{\Rightarrow }\mathit{x}{~}_{\mathit{\alpha }} \mathit{z} \underset{\square }{\iff }\exists {\left\{{\mathit{y}}_k\right\}}_{k=1,...,N}\in A / \left\{\begin{array}{c}d(\mathit{x},{\mathit{y}}_1)\le \alpha \\ d({\mathit{y}}_k,{\mathit{y}}_{k-1})\le \alpha \\ d({\mathit{y}}_N,\mathit{z})\le \alpha \end{array}\right., k=1,...,N-1$$

(2)

Given two elements C₁, C₂ ∈ M_2xN (Z), the distance $d_c$ between them is defined as:

$$\begin{gathered} d_c : M_{2\times N}\left(\mathbf{Z}\right)\times M_{2\times N}\left(\mathbf{Z}\right)\to \mathit{R} \\ { C}_1,C_2 \to d_c(C_1,C_2)=min(\left|x_1-y_1\right|+\left|x_2-y_2\right| \end{gathered}$$

(3)

where ($x_1,x_2$)∈ C₁ and ($y_1,y_2$) ∈ C_2.

3.2. Algorithm description

The main flowchart of the following algorithm is described in the Figure 4. It consists of three phases, which are explained below:

1.

Input elements:

• Pixel matrix set: Each matrix corresponds to the image associated with a cross-section of the pulp chamber, in ascending order of height. It is understood that the first image corresponds to the first plane of the pulp chamber that does not include the molar canals.
• Pixel dimension: For each set of images, the unit of distance is the measurement of the pixel side. In turn, the inter-planar distance corresponds to the height of the pixel. For the sets considered, we work with cubic pixels (voxel).
• First surface point: To allow the detection of the pixels of the pulp chamber surface in the first matrix, it is necessary to consider a pixel belonging to it. If this pixel is not provided by the user, a pixel is taken automatically in the manner indicated below.
• Control parameter k (threshold): In the images considered, there is a color difference between pixels that belong to the bony surface of the tooth and those that do not. The control parameter k makes it possible to establish a prior separation between the elements of the pixel matrix that potentially belong to the pulp chamber and those that do not.

2. 

Determination of the surface matrix:

Starting from the pixel matrix I of each plane, a new matrix M (surface matrix) is constructed such that:

$$M_{ij}= \left\{\begin{array}{c}1 if I_{ij} \ge k\\ 0 if I_{ij}<k\end{array}\right.$$

(4)

3. 

Calculation of pixel surfaces:

Given a non-zero element a of matrix M of coordinates (a₁, a₂), we define the equivalence class with representative a as ${\{M}_{ij}$ ∈ M / (i, j) ${~}_2$ ($a_1,a_2$) ∧ $M_{ij}\ne 0\}$. Each of the equivalence classes is called a surface.

In the first step of the algorithm a pixel of the pulp chamber is required. The first image, which is being worked with in the first step, is pointed out by an expert and is the closest to the molar canals. To pick up the pixel belonging to the pulp chamber automatically, a pixel from the center of the image is selected first. If the pixel value is greater than the threshold value, it is considered to belong to the pulp chamber and the algorithm is started. If the selected pixel does not exceed the threshold value, the next closest pixel to the center of the image, not previously considered, is taken and the process is repeated.

Since the first image is the one closest to the molar canals, we are under the assumption that all pixels in the pulp chamber, in that plane, belong to the same equivalence class.

To separate each of the equivalence classes, we start by taking a non-zero element in M (pixel), of coordinates P, from which a set S is defined, and this element is equated to zero in M. Subsequently, it is checked if exits the elements in M whose coordinates P′ fulfill

d(P,P′) ≤ 2 and verify that are not null. If so, they are introduced into S and the process is repeated for these new elements.

Once no more elements are added to S, the surface is considered complete and a new non-zero element is taken in M. The algorithm ends when there are no more non-zero elements in M to add to an equivalence class.

At the end of the process, we have a separation of M into its different surfaces. The flowchart for this phase is shown in Figure 5.

4. 

Identification of the pulp chamber.

If it is the first image, only that surface containing the first surface point given as input parameter is retained. If it is a subsequent image, surfaces far from those obtained in the previous iteration are eliminated. If $S_j^i$ denotes the j-th surface in the i image, then:

$$S_j^i\in pulp chamber \underset{\square }{\iff }\exists S_k^{i-1} / d_c(S_j^i,S_k^{i-1})\le 2$$

(5)

At the end of each iteration, if new pixels associated with the pulp chamber have been obtained and if there are still images of the input set to be evaluated, we proceed to work with the next matrix. Otherwise, the volume calculation is performed.

3.3. Volume calculation

The surface area of each pixel is the square of its side. Therefore, the value of the pulp

chamber image area corresponding to each plane is the total number of pixels a signed to the pulp chamber multiplied by the individual surface area of the pixels.

For the calculation of the pulp chamber volume, two approaches are considered. In one case, the areas of the lower planes are taken and multiplied by the value of the distance between images (planes). In the other, the areas of the upper planes are multiplied by the distance between images.

Thus, considering a total of N images, in each case N - 1 volume values are obtained, one for each interplanar distance, the corresponding approximation being the sum of the calculated volumes. In this way, the two approximations of the volume are the extremes of an interval whose central value is the value of the volume considered, and the radius is the uncertainty in the measurement.

To attempt to estimate the algorithm's complexity, a study of its runtime was conducted under average conditions. For this purpose, the set of images associated with one of the samples is taken as representative and resized, while maintaining their aspect ratio, to sizes two, three, four, and five times the original size. To study the average runtime, one hundred measurements are taken for each group of images corresponding to a size. Through least-squares adjustment, cubic growth is obtained (Figure 6).

Considering that the images processed in this case have a size smaller than 100x100 pixels, the average execution times obtained are just a few seconds. Thus, the functionality of the algorithm is proven as it does not require waiting times to verify if the result is satisfactory. On the other hand, the observed behavior with the increase in the number of pixels allows understanding the scalability of the algorithm to images with higher resolutions.

The average times were measured on a LAPTOP -AA7B78L5 Intel (R) Core (TM) i7-7500U CPU @ 2.70 GHz 2.90 GHz, with 8 GB RAM.

4. Results

The results with each method used in the study were obtained independently and by different analyzers and were supervised by two experts with extensive clinical experience to ensure the effectiveness of the procedure. The volumes of the pulp chamber according to the different methods used in the study are in Table 1. The data obtained with micro-CT were established as ideal and "true volumes", because to the high resolution of the images that were acquired due to the high dose of radiation used; those results were the most reliable.

The values obtained for samples 1, 2 and 4 were quite similar (Figure 2), while those for samples 3 and 5 were almost double (Figure 3). The range of the values was completely suitable and normal and was independent of the maxilla and place. From the CB-CT images manipulated by the specialist (CB-CT set threshold) the values were like those obtained by micro-CT except for samples 4 and 5 which were considerably lower. However, the volumes obtained with CB-CT automatic threshold were always much lower than those of micro-CT and CB-CT set threshold, except for sample 3. The volumes generated by the mathematical formula resulted reliable in 4 of the 5 samples analyzed, with an error of 3.65 % on average, and matched those obtained with micro-CT. Only in the sample number 5 the results calculated by algorithm clearly differed from those obtained with micro-CT.

From the point of view of the algorithm, it can only be said that the method of detection of the pulp chamber is performed automatically using as data the distance between planes, the initial plane, and the control parameter k. For the algorithm to work, it is also necessary to know the location of a pixel belonging to the pulp chamber of the first image, which can be given by the user or, otherwise, the central value of the image is taken. Based on this data, the pulp chamber is automatically detected, in the planes where it exists, and the corresponding volume. Therefore, from the perspective of the algorithm, there are no differences between sample 5 and the others (Figure 3), so it does not allow to justify the discrepancy between the volume provided by the micro-CT and the calculated.

Several authors have calculated the pulp chamber volume for different reasons and using different methods: measuring the volume of the pulp cavity from cone beam computed tomography (CBCT) images by comparing them with the micro-CT reference standard, using the results of different radiological units and comparing the 3D reconstructions using statistical methods ([13]); understanding internal anatomy relationships before undertaking endodontic therapy using artificial intelligence techniques by training the U-Net network for automatic pulp segmentation ([14,15]); for the estimation of human age by statistical methods ([16,17]). The results obtained by the segmentation algorithm developed in this work, from the axial CT images, are comparable to those obtained in other works and easily reproducible, the only condition being that the first plane containing parts of the pulp cavity must be pointed out by an expert.