Enhancement of Cone-Beam Computed Tomography Dental-Maxillofacial Images by Sampling Kantorovich Algorithm

In this paper, we establish a procedure for the enhancement of cone-beam computed tomography (CBCT) dental-maxillofacial images; this can be useful in order to face the problem of rapid prototyping, i.e., to generate a 3D printable file of a dental prosthesis. In the proposed procedure, a crucial role is played by the so-called sampling Kantorovich (SK) algorithm for the reconstruction and image noise reduction. For the latter algorithm, it has already been shown to be effective in the reconstruction and enhancement of real-world images affected by noise in connection to engineering and biomedical problems. The SK algorithm is given by an optimized implementation of the well-known sampling Kantorovich operators and their approximation properties. A comparison between CBTC images processed by the SK algorithm and other well-known methods of digital image processing known in the literature is also given. We finally remark that the above-treated topic has a strong multidisciplinary nature and involves concrete biomedical applications of mathematics. In this type of research, theoretical and experimental disciplines merge in order to find solutions to real-world problems.


Introduction
Here, we state some new applications of a certain family of approximation operators, known under sampling Kantorovich (SK) operators, which have been deeply studied in recent years in connection to both theoretical and practical aspects.
The SK operators have been introduced in [1] in a one-dimensional setting, and subsequently, they have been extended to the multivariate frame in [2]. The SK operators (based upon non-negative kernels) belong to the wide family of linear (positive) operators, which is one of the most studied topics in approximation theory (see, e.g., [3][4][5][6][7]).
Further, the SK operators also represent the L 1 -version of the well-known generalized sampling series introduced by P.L. Butzer and their school at Aachen in the 1980s (see, e.g., [8,9]).
The new application that we study is related to the reconstruction and the enhancement of dental cone-beam computed tomography (CBCT) images (or cone-beam computed tomography, [10][11][12]), in connection to the problem of prototyping. In this regard, a crucial role is played by the SK operators and their implementation for image processing.
Dental cone-beam computed tomography (CBCT) is a very recent technology that allows for the production of maxillofacial non-invasive images, finding applications in the problem of rapid prototyping, i.e., prosthesis prototyping by means of a 3D printer.
In particular, to support the implantation of the dental prosthesis, it is helpful to achieve a 3D-printed model of the jaw of the region of interest (ROI), i.e., the area where the surgeon is going to perform the dental implant.
In the procedure that allows passing from CBTC images to the 3D prototyping (which works on bidimensional images), the application of the SK operators to the original CBCT image is the first step; this allows the image noise to be reduced (since this algorithm works as a low-pass filter), consequently improving the quality of the reconstructed image. Subsequently, a 3D reconstruction model for the generation of 3D printable files is also applied. In the first step of the above procedure, we can also consider other reconstruction methods known in the literature (bilinear B-spline and bicubic B-spline) in place of the SK algorithm, in order to establish what is that one that produces the best reconstruction of the CBCT images for this specific problem.
The comparison is made using the well-known peak signal-to-noise ratio (PSNR), an index commonly used to evaluate the quality and similarity of digitally processed images.
The present work suggests a solution to the problem of rapid prototyping for the automatic generation of 3D dental prosthesis. In the above research, it revealed to be fundamental the use of CBCT images provided by the second author, for the numerical test of the proposed procedure. To the best of the authors' knowledge, this is the first attempt to apply techniques of digital image processing based on sampling-type operators to the above, very current, problem, and this represents the main novelty of this research work. Then, this can be certainly classified as a "new direction" of "approximation theory", in which the computational aspects of a mathematical theory play a crucial role and merge in order to find solutions to real-world problems in the biomedical field. On the other hand, the idea of using an optimized implementation of certain theoretical (mathematical or, more precisely, approximation) results in order to establish concrete applications that suitably fit with current research trends that follow innovative way strictly related to the "problem solving" and to the development of new technologies. In this sense, the experimental results established here follow the latter mentioned direction.
The present research can also open the way to several future works, such as possible quantitative analysis of the application here proposed, and also several comparisons with other well-known noise filters and methods of image processing in order to establish which of these can produce the best results for 3D prototyping.

Materials and Methods
Below, a detailed description of the materials and methods used in the present research is presented.

Cone-Beam Computed Tomography Dental-Maxillofacial Images
Dental cone-beam computed tomography (CBCT) technology emerged in 1995 when Italian inventors A. Tacconi and P. Mozzo introduced the first maxillofacial imaging device, namely, the NewTom DVT 9000. The latter is a diagnostic device based on X-rays. X-ray images (or radiograms) depend on the attenuation of the X-ray beam that crosses the body; if the density of the structure increases, then the X-rays intensity becomes attenuated, and vice versa.
Dense objects produce an opacity on the radiogram (hence, they are called radioopaque), while the objects having low density (called radio-lucent) appear as dark areas in radiograms.
CBCT is a special kind of CT (see, e.g., Figure 1) whose geometry and size vary according to the position (standing, sitting, or supine position) that the patients assume during the diagnostic examination. Furthermore, the size of the field of view (FOV) is variable and depends on the size and shape of the detector, the beam projection geometry, and the ability to collimate the beam. In CBCT, the FOV is the most important scanning parameter affecting the patient radiation dose and image quality, and therefore, it should be selected according to clinical needs. CBCT uses an X-ray source that turns in a unique, complete scan (360 degrees) around the investigated object such that the X-ray beam (that can be continuous or pulsed) has a cone shape. In this way, in a single scan, it is possible to acquire a series of two-dimensional images of the anatomical region of interest (ROI).
The number of projections collected during the rotation depends on the detector frame rate and usually varies from 128 to 1024, giving different bases for the image reconstruction. Image quality is affected by the number of projection images, an increased number of samples basically resulting in higher image quality. The X-ray beam emitted by CBCT crosses the patient and hits a bidimensional detector placed on the opposite side of the object relative to the X-ray source that detects the absorption data. A process of data analysis and reconstruction follows to obtain a series of images diagnostically valid in every observation plane.
Image reconstruction in CBCT is a mathematical process in which the measured projection raw data allow the representation of the anatomy of the patient by means of images. The attenuation value distribution detected from the patient is presented in the 3D CBCT image data by a grayscale. Indeed, the density of tissues is expressed by grayscale voxels built on the basis of the number of Hounsfield Units (HUs). The Hounsfield number is a value proportional to the absorption of the X-rays and connected to the density of the tissue.
We emphasize that the main importance of the cone-beam images is the reconstruction of a 3D image of anatomical structures with lower radiological invasiveness than the classic CT. Indeed, we can assume that an average radiation dose for a cone-beam CT of the jaws taken for implant purposes is approximately 130 µSv, where Sv stands for Sievert; it denotes a derived unit of ionizing radiation dose in the international system of unit (SI), and it is a measure of the health effect of low levels of ionizing radiation on the human body. This means that it is much less invasive than a CT (a CT of the head is approximately 860 µSv) and that the amount of radiation received for a cone-beam CT is minimal in comparison to what we receive being alive on the earth (the level is similar to the amount received while being alive on the earth for approximately 16 days).
Moreover, thanks to CBCT, the tomography slice can be displayed, at the same time, along a number of different directions as, for instance, the three orthogonal planes, i.e., the axial, the sagittal, and the coronal planes.
The properties of CBCT are well suited for dental examinations in the maxillofacial region where high spatial resolution of hard tissue is, typically, of primary interest. As for other imaging techniques, cone-beam images are also typically affected by artifacts, in this specific case, due to the hardening of the beam and to the partial volume effect (PVE). The first one is caused by the absorption of the low energy radiation by the tissue with elevated atomic number (as bone), while the second occurs when there are structures in a voxel with very different densities.
In the image, the different density of the achieved voxel is represented by a mean that does not correspond to any real density [13].
The presence of artifacts and noise, justify the applications of postprocessing methods to CBCT images, such as enhancement algorithms and low-pass (noise) filters.

The Sampling Kantorovich Operators
In this section, we recall the main definitions related to the theory of sampling Kantorovich (SK) operators, which are on the basis of the SK algorithm.
Here, we denote by t k = t k 1 , ..., t k n a vector in which each for every i = 1, ..., n. Note that, the elements of (t k i ) k i ∈Z are not necessary equally spaced. The symbol Π n denotes the sequence (t k ) k∈Z n taken from the whole (symmetric) R n .
A function χ : R n → R is said to be a kernel if it satisfies the following properties: where · 2 denotes the usual Euclidean norm. We now recall the definition of the multivariate sampling Kantorovich operators introduced in [2]. Define where f : R n → R is a locally integrable function such that the above series is convergent for every x ∈ R n , where The operators in (1) have been originally introduced in [1] in the univariate setting as an L p -extension of the well-known generalized sampling series. Below, we recall the main approximation results which justify the applications discussed in the next sections. We have the following: where C 0 (R n ) is the space of all continuous and bounded functions f : R n → R.
In particular, if f ∈ C(R n ), i.e., f belongs to the space of the uniformly continuous and bounded function, then lim where · ∞ denotes the usual sup-norm.

Theorem 2.
For every f ∈ L p (R n ), 1 ≤ p < +∞, the following inequality holds: where · p denotes the usual p-norm. Moreover, we also have Note that Theorem 2 is a direct consequence of a modular convergence result established in the general context of the Orlicz spaces ([2]).
The above operators, other than the above convergence results, have been studied in depth under several theoretical aspects. For instance, the qualitative and quantitative order of approximation (with respect to different norms) has been established in [14].
In order to understand what their best performance can be, saturation and inverse theorems of approximation have been established in [15,16]. In particular, in [16], the saturation order with respect to the L 1 -norm has been established in a subspace of L 1 , resorting to the so-called Fourier transform method.
A crucial fact in the applications of the sampling Kantorovich operators that will be discussed in the next sections resides in the choice of the kernel functions. About this problem, a detailed discussion, together with numerical tests, can be found in [17].
Multivariate kernel functions can be considered in different ways. For instance, one can take into consideration radial kernels, i.e., functions for which the value depends on the Euclidean norm of the argument only. An example of such a kernel can be given, for example, by the Bochner-Riesz kernel, defined as follows: for x ∈ R n , where α > (n − 1)/2, B λ is the Bessel function of order λ and Γ is the Euler function. For more details about this matter, see, e.g., [8].
Further, multivariate kernels can also be constructed using the product of univariate kernels, see, e.g., [2,8]. We briefly describe this procedure in the case of sampling Kantorovich operators in which we consider a uniform sampling scheme, i.e., the sequence t k = k, k ∈ Z n . Now, we denote by χ 1 , ..., χ n , the univariate functions χ i : R → R, χ i ∈ L 1 (R), satisfying the following assumptions: for some β > 0, χ i is bounded in a neighborhood of the origin and for every u ∈ R, for i = 1, ..., n. Now, setting χ(x) := ∏ n i=1 χ i (x i ), x = (x 1 , ..., x n ) ∈ R n , it is easy to prove that χ is a multivariate kernel for the operators K w satisfying all the assumptions of our theory, again see, e.g., [8].
As a first example, consider the univariate Fejér kernel defined by F(x) := 1 2 sinc 2 x 2 , x ∈ R, where the sinc function is given by The Fejér kernel is a symmetric (more precisely even) kernel. Then, we can define ., x n ) ∈ R n , the multivariate Fejér kernel satisfying the condition upon a multivariate kernel. The function F n is an example of a kernel with unbounded support; then, to evaluate our sampling Kantorovich series at any given x ∈ R n , we need an infinite number of mean values w n R w k f (u) du. This problem can be overcome if we consider function f having compact support, while, in the case of a function having unbounded support, the infinite sampling series must be truncated to a finite one, which leads to the so-called truncation error. In order to avoid the latter approximation error, one can take duration limited kernels, i.e., functions χ with compact support. Useful examples of kernels with compact support can be constructed using the well-known central B-spline of order k ∈ N ( [18]), defined by where the function (x) + := max{x, 0} denotes the positive part of x ∈ R (see [1] again). Furthermore, the central B-spline of order k are examples of one-dimensional symmetric kernels. Hence, by M n k (x) := ∏ n i=1 M k (x i ), x = (x 1 , ..., x n ) ∈ R n , we denote the multivariate B-spline kernel of order k ∈ N.
Finally, other important examples of univariate kernels are given by the Jacksontype kernels, defined by J k (x) = c k sinc 2k x 2kπα , x ∈ R, with k ∈ N, α ≥ 1, where the normalization coefficients c k are given by In a similar manner, using the previous procedure, we can construct a multivariate version of Jackson-type kernels. For more details about Jackson-type kernels, and for other useful (not necessarily symmetric) examples of kernels, see, e.g., [19][20][21][22][23][24][25][26].

Enhancement of CBCT Dental-Maxillofacial Images by the SK Algorithm
In this section, we show applications of the bivariate sampling Kantorovich operators to image reconstruction (see, e.g., [27]), with particular attention to the problem introduced in Section 2.1, i.e., the enhancement of dental-maxillofacial CBCT images.
The above function I(x, y) is defined in such a way that, to every pixel (i, j), it is associated with the corresponding gray level a ij .
We can now consider the approximation of the original image by the bivariate sampling Kantorovich operators (K w I) w>0 based upon some two-dimensional kernel χ.
Then, in order to obtain a new image (matrix) that approximates pointwise and in L p -sense the original one, it is sufficient to sample K w I (for some w > 0) with a fixed sampling rate. In particular, we can reconstruct the approximating images (matrices) taking into consideration different sampling rates, and this is possible since we know the analytic expression of K w I.
A pseudo-code related to an optimized implementation of the SK algorithm can be found in [28].
The analysis of the reconstructions obtained by the application of the SK algorithm to CBCT images will be performed using the well-known peak signal-to-noise ratio (PSNR) index, defined as follows: where f max is the maximum possible value of the signal, or function f (the full-scale value), and MSE is the standard mean square error, defined in the (bounded) support D ⊂ R n of f , as follows: where f is the original signal, and f r being the reconstructed version of f . The PSNR is extensively used in image analysis and processing to evaluate, for example, the rate of similarity of two images after a reconstruction process, see, e.g., [29]. Of course, in the case of digital images, we use a discrete version of both MSE and PSNR. Note that the lower the MSE is, the more the f r reconstructs qualitatively better the original image, which implies higher values of the PSNR are representative of a better reconstruction.

Numerical Results
Numerical simulations were performed to evaluate the quality of reconstructions using different methods. In particular, we compared the results achieved by bilinear and bicubic algorithms with the ones coming from sampling Kantorovich operators. We used bilinear and bicubic interpolation methods, as defined in [17], i.e., bilinear and bicubic B-splines. All the simulations were operated in MATLAB©. Two CBCT sets of images were used: the first relative to an entire arch of the upper jaw of a patient, and the second relative to an entire arch of the lower jaw of another individual. The code used for the reconstruction was entirely written in MATLAB©. All the images used are in 16 bits codification, and no compression on the data was performed, exploiting the whole dynamical range of the acquisitions and avoiding any windowing procedure. The SK algorithm was applied with the parameter w = 20, the scale factor r = 2, and using the Jackson-type kernel with k = 3. The choice of the kernel and of parameters for the implementation of the SK algorithm were made according to the results established in [17]. Example images extracted (at different depths) from the lower and upper dental arches are shown in Figures 2 and 3, respectively.  From a visual point of view, all the images appear really similar to each other; however, there are some differences that are not visible to the naked eye. In Figure 4, the differences between the SK operator reconstructions and the bilinear and bicubic are shown. Most of the differences are localized along the borders of the dental arch, a particularly meaningful aspect for the reconstruction of the correct shape of the dental implant. The differences are mostly concentrated in the dental profile, an aspect that is particularly important for a correct implant design.
Moreover, the numerical estimation of similarity by means of PSNR is reported in Table 1, where the original CBCT images are taken as reference. The calculation of PSNR was performed on axial images using the native psnr() MATLAB function, based on the definition given in (4). The achieved values are not operator dependent. Table 1. Numerical evaluation of PSNR for bilinear B-spline, bicubic B-spline, and SK operators reconstructions. The first column indicates the number of sections considered in the whole scan. The results for other images of the scan are qualitatively the same as the ones shown in the table, with SK performing anyway better than the other methods (SK used parameters are w = 20, r = 2, k = 3, where k is the parameter of the Jackson kernel).

Discussion and Conclusions
From the numerical results of Table 1, we can deduce that the SK reconstructions offer better results because they have a higher value of the PSNR than the other considered methods.
To extend the analysis concerning the stability of the reconstruction even in the case of noisy images, another scan was considered. Achieved results are shown in Figure 5. Moreover, the numerical estimation of similarity between noisy images and their reconstructions by means of PSNR is reported in Table 2, where the original CBCT image is taken as reference. We highlight again that the proposed method seems to be performing even in the case of original images affected by noise and with very low quality.

Future Developments: Applications of 3D Prototyping
In this section, we show the procedure that allows CBCT images to be printed using 3D printing technique. In particular, we want to print the region of interest (ROI), i.e., the area where the surgeon is going to perform the implant, see Figure 6.
A dental implant is a surgical component that interfaces with the bone of the jaw or skull to support a dental prosthesis. The basis for modern dental implants is a biologic process called osseointegration, which consists of the direct structural and functional connection between living bone and the surface of a loadbearing artificial implant. This is possible thanks to some materials, as titanium, that can form an intimate bond to the bone. The implant texture is first placed so that it is likely to osseointegrate. A variable amount of healing time is required for osseointegration before either the dental prosthetic (a tooth, bridge, or denture) is attached to the implant. The standard process to realize an implant is quite annoying for the patient, for the fact that the realization of a dental mold of the mouth is needed. The dental mold is necessary both for the design of the implant as for the manufacturing of the surgical guide that will be placed between the two teeth inside the ROI. This guide is used to hold the drill during the gingival trepanation that is necessary to put a pin on which screwing an artificial tooth. The goal of our procedure is to get the prototype of the dental mold directly by CBCT images using the 3D printer; as a consequence, as above said, it is functional for the design of the surgical guide.
In this way, there is no need to realize the dental mold, while the dose of the radiation administered to the patient is the same since CBCT exam is a standard clinical preoperational exam.
Following these results, an example of the future development of the applications considered in the present paper will be the realization of the prototype realized by 3D cone-beam images, after the improvement achieved by the application of the SK algorithm. Some preliminary applications in this sense are given in Figures 7 and 8, where 3D models are given under various visual points of view.
This can be useful since, as previously discussed, in the images processed by the SK algorithm, the noise is reduced, and the quality of the 3D prototype seems to be better, both from an anatomical point of view and because the 3D printing process needs a sufficiently smooth surface, a condition that is guaranteed by the smoothing procedure inherent to the SK algorithm.  Finally, in order to support the above claim, a quantitative analysis should be performed. The main problem arising in such a numerical test is due to the application of a suitable registration procedure for 3D images, which is a fairly complex study that could be carried out in the future.
Author Contributions: The authors D.C., P.P., M.S. and G.V. contributed equally to this work (concept, writing, methodology, review and editing). All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: