Spectral photon counting for panoramic dental imaging

Panoramic x-ray imaging is a versatile, low-dose imaging tool, which is routinely used for dental applications. In this work, we explore a further improvement of the concept by introducing recently developed spectral photon-counting detector technology into a conventional panoramic imaging unit. In addition we adapt spectral material decomposition algorithms to panoramic imaging needs. Finally, we provide first experimental results, demonstrating decomposition of an anthropomorphic head phantom into soft tissue and dentin basis material panoramic images, while keeping the noise level acceptable using regularization approaches. The obtained results reveal a potential benefit of spectral photon-counting technology also for dental imaging applications.


Introduction
Panoramic is a useful, low-dose imaging modality for dental radiography (Gijbels et al 2005). One way to generate a panoramic image is by summing partially overlapping frames together (McDavid et al 1995). The shift from one frame to the next is determined by the movement and rotation of the imaging system. Panoramic imaging provides an overall view of the maxilla, the mandible, and the teeth while being able to show many lesions in this area (Ianucci and Howerton 2012). Panoramic imaging has also been used for detecting caries (Kamburoglu et al 2012), periodontitis (Nardi et al 2018), and maxillary sinusitis (Ohashi et al 2016). However, for all these conditions the intra or extra observer agreement in diagnosis is not optimal. One reason for sub-optimal diagnostic accuracy is the superimposition of soft and hard tissues in this imaging modality. This could eventually be overcome by recently developed spectral x-ray imaging approaches (Mechlem et al 2018, Flohr et al 2020, Shi et al 2020. These approaches have demonstrated superior diagnostic capability by employing material decomposition approaches based on spectral detection and the ability to detect features that would remain hidden in a conventional image. One challenge, however, is that for panoramic imaging, the photon statistics for an individual frame are very low. In other words, they have a very low number of average counts per frame and thus suffer from bias. This implies that the pre-existing spectral projection domain algorithms, e.g. (Roessl and Proksa 2007, Alvarez 2011, Mechlem et al 2018, cannot be directly applied to panoramic imaging data, as applying these algorithms to very low photon statistics observed in panoramic imaging would introduce statistical bias (Mechlem et al 2018). The decomposition is a nonlinear operation and the noise is propagated through the decomposition, leading to an amplified noise in the output of the decomposition .
One way to do spectral panoramic x-ray imaging is to calculate the ratio of bins during imaging and compare that ratio to one obtained from reference phantom scan , Langlais et al 2015. However, the output of this algorithm is one unitless number i.e. not a line-integral of two or more basis materials as e.g. by an algorithm used for spectral angiography (Mechlem et al 2018). Furthermore, this Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. algorithm does not thoroughly address the spectral responses of individual detector pixels.
In this paper, we present a first quantitative material decomposition algorithm for panoramic imaging. We show how to adapt an existing material decomposition algorithm to address the requirements of panoramic imaging. The results of the decomposition are two basis material panoramic images. The performance of this algorithm is demonstrated in an experimental setup using an anthropomorphic head phantom.

Material decomposition
The concept of using low-order polynomials to model line-integrals in spectral x-ray imaging was first demonstrated by Cardinal and Fenster (Cardinal and Fenster 1990). In their work, the decomposition for two spectral measurements and materials was done by having the basis material line-integrals (A i 0 , A i 1 ) for detector pixel i as a function of the low and high energy log-signals (l i l , l i h ) and parameters determined in This kind of formulation enables a fast and accurate decomposition in the case where there are the same number of basis materials and spectral measurements by simply evaluating these functions at the log-signals of the measurements. However, it is not possible to incorporate any kind of regularization for this closedform solution.
If a maximum likelihood (ML) estimation is desired, then the inverted relation (Mechlem et al 2018) can be used. Here, the expected number of counts y i ŝ for pixel i and spectral measurement s is defined as: where b i s is the flatfield signal for pixel i and spectral measurement s. The polynomial model P i s used in this study assumes two basis materials: where c (pixel and spectral measurement subscripts dropped for simplicity) are the fit coefficients. Initially, a least-squares fit to the calibration measurements is carried out to find the fit parameters c is : where l ik s are linearized calibration measurements, A ik are the line-integrals, and w k the weights with k being the calibration measurement index. The weights are assigned according to the standard deviation of the corresponding calibration measurement. During calibration, the line-integrals (i.e. thicknesses) of the materials are obtained by measuring these values from the physical calibration phantom. Following the forward model calibration (equation (4)), a maximum likelihood decomposition is carried out to find the combination of basis materials thicknesses A i maximizing the likelihood given the spectral measurements. Assuming that the spectral measurements follow independent Poisson statistics, it is computationally more efficient to minimize the corresponding negative log-likelihood: In equation (5), the expected number of counts are a function of the basis material line-integrals. This enables maximum likelihood estimation of the lineintegrals. Generating the basis material images (equation (5)) leads to noise amplification and a degradation of contrast to noise ratio. Both of these problems can be remedied by denoising which takes into account the anti-correlated noise in basis material images. The where -A ( ) is the sum of negative log-likelihoods over all pixels and spectral measurements. A ( ) is the regularization term which is responsible for the denoising. q denotes patch index. Dictionary D contains a large number of small, high-quality images of typical structures of relevant anatomy. These small images in a dictionary are called atoms. The image to be denoised is divided into small patches, which are still larger than atoms. The idea is then to find a sparse linear combination α q of atoms that represent the image patch E q A to be denoised. α q and A are each optimized in turn so that one is optimized while the other remains constant. j q is the whitening transformation matrix converting E q A to new bases so that the noise is uncorrelated across the spectral channels.

Panoramic imaging
In panoramic imaging the imaging system follows an elliptical trajectory around the patient so that the detector is on the front of the patient and the tube is on the back. The panoramic image is formed by summing partially overlapping frames together along this trajectory (Noujeim et al 2011). Figure 1 shows a working principle of panoramic imaging. The sharp layer typically aligns with the dental arch. Everything outside the sharp layer gets blurred while the details in the sharp layer remain visible. Additionally, it is important that one ray path crosses the dental arch only once (so that the teeth don't overlap) and that the path is roughly perpendicular to the arch. The process of generating the panoramic image is sometimes called reconstruction even though no backprojection is done, as opposed to computed tomography (Ianucci and Howerton 2012).
Post processing is typically applied after the panoramic image has been formed. Typically, the purpose of this post processing is to remove noise and adjust brightness and contrast so that more details are visible in single window level setting. (Gormez and Yilmaz 2009).

Experimental measurement 2.3.1. Phantom
An anthropomorphic head phantom (Erler-Zimmer GmbH, Lauf, Germany) was used in this study as seen in figure 2. The phantom features a realistically attenuating skull, a dental bridge and a dental implant all encapsulated in plastic. These kinds of phantoms are used in dental x-ray unit production quality assurance among other things.

Experimental setup
A DC-Vela photon counting detector (Varex Imaging Corporation, Salt Lake City, USA) was fitted to a Planmeca Promax Mid imaging unit (Planmeca Oy, Helsinki, Finland) as shown in figure 2. The detector has two adjustable energy thresholds, a charge-sharing correction, a 0.75 mm thick CdTe sensor layer and a maximum frame rate of 300 fps when using two thresholds. The dimensions of the sensor are 64 (width) pixels by 1541 (height) pixels and the detector has a pixel size of 100 μm. The adapter holding the sensor was designed so that the source-to-image distance (SID) and the imaging geometry is the same as in Planmeca's standard imaging units.
The panoramic reconstruction was carried out using Planmeca's proprietary back end software. The same back end software is used in Planmeca's commercially available units.

Acquisition parameters
Prior to the phantom measurement, 36 calibration measurements with different thicknesses of titanium and polyoxymethylene (POM) were performed. The thicknesses are listed in table 1. Titanium and POM were chosen for basis materials as their attenuation profiles (i.e. effective atomic number and density) (Cullen et al 1997) are close to those of dentin (Zenóbio et al 2011) and soft tissue (Hubbell and Seltzer 1995), respectively. The latter pair reflects tissue types typical to jaw and lower head.
Acquisition Parameters are shown in table 2. The dose-area product (DAP) was measured with Ker-maX Plus IDP 120-131 HS sensor (IBA Dosimetry GmbH, Schwarzenbruck, Germany). According to a 2019 diagnostic reference study (Hodolli et al 2019), the mean DAP of 21 panoramic imaging units for adult patients was 74.1 mGy(cm) 2 with a standard deviation of 28.9 mGy(cm) 2 . The DAP of 72.1 mGy(cm) 2 used in this study is very close to the mean and therefore clinically relevant. The charge-sharing correction was on for all measurements. Different tube currents were used for calibration and imaging (c.f. table 1). During calibration the detector is exposed to flatfield radiation and a higher current would induce pulse pileup in the detector application-specific integrated circuit (ASIC). Thresholds were chosen so that the low energy threshold (15 keV) captures all counts and high threshold (36 keV) has roughly half of the counts of the low threshold for the most attenuating calibration step. This was not an explicit optimization for any specific task as in (Wang and Pelc 2011) but it was manually found to be a working solution. The number of photon counts during calibration was of similar magnitude in this study (table 2) and a previous study (Mechlem et al 2018). However, in this study the exposure time was not increased for thick calibration measurements. The number of flatfield photon counts during imaging in this study is low (1.2 × 10 4 ). In a similar material decomposition approach (Mechlem et al 2018) statistical bias starts to quickly increase below 10 3 counts, which is well below the number of flatfield counts in our study. Additionally, the use of regularization has been demonstrated to lower the bias limit . However, the previous study (Mechlem et al 2018) used different photon spectra and a differently sized sample, which alters the number of detected counts and affects the noise properties and the bias limit.

Material decomposition workflow in panoramic imaging
The algorithm presented in 2.1 takes in 2D native detector frames. However, native detector frames cannot be used to directly decompose panoramic images as described in section 2.2 as the native frames have different dimensions and amount of photon counts in contrast to a panoramic image. Furthermore, the photon statistics in the native frames are low. Such low statistics would introduce a significant bias in the estimation (Mechlem et al 2018). Thus, some additional steps are required to be able to get quantitative material decomposition images. The workflow to generate spectral panoramic images is shown in figure 3. The first step was to acquire calibration frames in the native detector size 5000 native (64 × 1541) frames were acquired for each of the 25 calibration measurements. The mean of each calibration measurement is taken to suppress Poisson noise. The panoramic reconstruction in this work was done by summing 3700 partially overlapping frames. To generate a virtual panoramic image, instead of Figure 3. Spectral panoramic imaging workflow. In virtual panoramic images the detector frames are summed in similar fashion to actual panoramic. Thus, the value of each pixel in the output image of a virtual panoramic is the sum from the same pixels in the same frames as in actual panoramic i.e. both virtual and actual panoramic images have the same pixelwise spectral response. This ensures that the model works correctly for actual panoramic images. Calibration needs to be run only once and the same coefficients can be used for different panoramic measurements. using partially overlapping individual frames, the mean frame was used as every frame. A virtual panoramic image was generated for each calibration measurement. The size of these virtual panoramic images was 3390 × 1541, matching the size of actual panoramic images. The calibration was then carried out for these images by fitting the polynomial coefficients according to equation (4). Panoramic images are acquired normally. These images are decomposed as any 2D image, using equation (5). Regularized decomposition is done using equation (6) as explained in section 2.1. Finally, virtual monoenergetic images are generated by multiplying the decomposition results and their respective linear attenuation coefficients and summing together these two multiplications.

Results
The accuracy of the calibration coefficients is assessed by using them to decompose a set of test points having different thicknesses of POM and titanium than those used in the calibration (c.f. table 1). The calibration coefficients were first calculated using equation (4). The test points were then decomposed using equation (5). The mean of each test point was taken from over the whole virtual panoramic image, excluding some bad pixels, mostly near the edges of the image. The true values versus the results of the test decomposition are shown in figure 4. Some of the decomposition points (>0.6 cm of titanium, >14 cm of POM) in figure 4 actually exceed the maximum thicknesses of the calibration data points in table 1. Thus, it is not feasible to expect the model to yield correct values here. For the last calibration point (14 cm POM, 6.1 mm Ti) the radiation had attenuated by a factor of 4.1 × 10 −4 for bin low and 1.7 × 10 −3 for bin high. During imaging the total number of counts expected for flatfield is (c.f. table 2) 1.2 × 10 4 . Thus, only a couple of counts are expected with such attenuation. Having an attenuation which yields very small number of counts will be very noisy and suffer from rounding error due to the discrete nature of the photons. Note that this does not affect calibration test decomposition in figure 4, as there are more flatfield counts during calibration (c.f. table 2) than during imaging.
A 2D histogram of the output of equation (5) applied to the panoramic image obtained from the anthropomorphic head phantom is shown in figure 5. The histogram has 101 bins in each direction. A Gaussian filtering with a sigma of 2 was applied to the histogram in order to the make the prominent features clearer. The test calibration points with indices from 7-9, 13-15 and 19-21 (c.f. Figure 4) span the region of the line-integral values the decomposition of the head phantom yields. In other words, the accuracy of these 9 points is relevant for the sample studied and the reliability of the decomposition of the head phantom. The maximum absolute error for this set of test decomposition points is 0.028 cm of titanium (index = 8), with the rest of the errors being around 0.01 cm in figure 4. For POM the maximum absolute error for this set is 0.55 cm (index = 21), with the rest of the errors being less than 0.4 cm.
Subsequently, a change of basis from POM and titanium was performed into dentin and soft tissue using coefficients from the EPDL97 database (Cullen et al 1997). The output of the algorithm described in section 2.4 is shown in 6, where the dentin (b) and the soft tissue (c) images in figure 6 have units of length (cm). Figure 6 also contains the flatfield corrected and linearized all counts image (a) as a reference and as an indication of the total attenuation. Figure 7 contains the 50 keV virtual monoenergetic image which is unitless. Two regions of interest (ROIs) were chosen from the monoenergetic image in figure 7 so that ROI 1 had considerable hard tissue component while ROI 2 had only soft tissue component. This is reflected on plots in figure 7. More specifically, ROI 1 has considerably larger absolute values in attenuation than ROI 2 and the ratio of the absolute values decreases as the function of energy.
The soft tissue component in figure 6 is very homogeneous. The reason for this is that the phantom is a skull fully encapsulated in plastic. The skull displaces some plastic so that larger values in the dentin image (b) correspond to locally smaller values in the soft tissue image (c), compare e.g. the lower teeth in figures 6(b) and (c) below the red dashed ellipse. In a living human patient there would be air cavities e.g. in the mouth, the nose and the throat. The increased soft tissue values around the center of the image around the horizontal direction are a result of actually having more plastic in the beam path in these directions i.e., the increase in the values in this area is not an artefact. This increased thickness is visible in figure 7 and in the soft tissue image in 6, but not in the dentin image in figure 6.
A threshold was applied for inpainting very low photon statistics in the panoramic image before decomposition. Here, a threshold of 10 counts /bin/ pixel was applied to the decomposition output as a very low number of counts is typically caused by a metal object in the beam path. Thick metals are not decomposed correctly when there is also a soft tissue or POM component present, e.g. test calibration points 23, 24, and 28-30 in figure 4. Thus, the dental bridge over the bottom left molars and the implant on a bottom right molar in figure 6 appear almost monotonously white in the dentin image (b) and black in the soft tissue image (c).
The white blob within the red dashed ellipse in the linearized all counts image (a) figure 6 is most likely caused by an overlap of the mandible from the other side with the teeth or the maxilla. Thus, for x-ray paths in this area the attenuation is increased. However, this area is actually darker in figure 6 in the dentin image (b) and again brighter in the soft tissue image (c). This is an indication that the estimation might not work correctly for this area as it does not make sense that the dentin component in this area is smaller than in the surrounding area. In other words, the situation is somewhat similar as in the case of metals, but the number of counts in this area is still above the threshold.
The skull and the teeth are visible and have positive values in the dentin image in figure 6, while the values between the spine and the jaw are close to zero (e.g. ROI 2 in figure 7), which is to be expected as there is no bone along those x-ray paths. Most details, e.g. the pulps of the teeth, are still visible in the dentin component in figure 6 similarly to the all counts image in the same figure. However, the contrast for some of the very finest details has been reduced in the dentin image in figure 6 (b) in comparison to the linearized all counts image (c.f. the blue ellipses in (a) and (b) in figure 6).
Areas within the orange and green dashed circles have almost the same intensity the dentin image (b) in figure 6 but in the linearized all counts image (a) the area within the green dashed circle is darker than the within the orange circle. The latter trend is the same for the soft tissue image (c). In other words, the mandible seems to have somewhat uniform thickness in the dentin component while there is less soft tissue component (plastic in the phantom) along the line-integral paths within the green dashed circle than within the orange circle, according to the output of the algorithm. Recall that the flatfield corrected image represents the total attenuation i.e., both the dentin and the soft tissue components contribute to it. Thus, it makes sense that the mandible area in the linearized all counts image figure 6(a) does not have uniform intensity. The mandible area is mentioned as an example here, but the trend is visible over the whole image.
As the decomposition is carried out for the reconstructed panoramic image, some anatomic noise is to be expected. In other words, several frames having slightly different paths through the anatomy contribute to one column in the final panoramic image. Most of the bone (the mandible, the maxilla, and the teeth) are withing the sharp layer and do not need to be considered. However, the spine is not within the sharp layer. Furthermore, sometimes the opposing mandible overlaps with the anatomy in the sharp layer and might cause an error. Finally, most of the soft tissue is expected to be outside of the sharp layer, and needs to be considered as well.
The width of the detector (64 pixels) would be an intuitive window size when measuring the horizontal change of thickness in a panoramic image. However, in reality the majority of the signal is carried in the center pixels of the frame, with the edge pixels having only a minor contribution. This is the case due to the collimator blades blocking the majority of the signal at the edges. Thus, we used the full width at half maximum of the mean flatfield frame, 30 pixels, as the width to estimate anatomical noise.
The anatomical noise of the soft tissue component was quantified by estimating the change of thickness in the soft tissue image in figure 6. The general trend of the image along the horizontal dimension is smooth and almost constant. The strongest change is along the center area, where the maximum thickness change is around 1 cm / 30 pix.
The dentin image in figure 6 was used to measure the change of the bone signal for spine and the opposing mandible. For the spine, the maximum change was around 0.1 cm / 30 pixels, which was measured from the narrow area below the mandible. For the opposing mandible, the change was 0.25 cm / 30 pixels for the more prominent right side of the image, and 0.05 cm / 30 pix for the left side. This was measured by plotting the edge between the teeth.

Discussion
The workflow described in section 2.4 enables the decomposition of panoramic images as 2D images at clinically relevant dose levels. However, the method does lose some spectral sensitivity by summing pixels along the horizontal dimension. Different pixels in a physical detector have slightly varying energy responses and summing them together widens spectral overlap across the energy bins.
Regarding anatomical noise, the center pixels for the soft tissue in figure 6 were identified as an area where the noise is the highest. For bone, the border with opposing mandible and the spine were identified as having some noise. The implication here is that the neighbouring columns have also contributed to the absolute value when looking at the soft tissue or dentin images at the respective areas. Finally, details outside the sharp layer will be blurred. However, this is the case for panoramic in general, not just the spectral algorithm presented in this paper.
Even after the calibration, the model presented in section 2.1 has finite accuracy in decomposition as indicated by figure 4. Errors get larger if the material thicknesses in imaging are larger than in the calibration. Furthermore, the model has not been designed to account for scatter. Even though the calibration materials were placed in contact with the collimator, due to the small source to image distance, some amount of scatter still reached the detector. For increasing thicknesses, more photons interact via Compton or Rayleigh scattering, while the number of unattenuated photons gets exponentially smaller. The model assumes that photons which interact with the object being imaged do not cause counts in the detector. This is strictly speaking never the case, but issues appear to arise only for the strongly attenuating test decomposition points in figure 4. In contrast, when observing the distribution of panoramic decomposition results in figure 5 we see that a vast majority of the points are <12 cm of POM, <0.2 cm titanium. As the test decomposition points in this regime are decomposed quite accurately in figure 4, a low decomposition error can also be assumed for the sample scan. Finally, it should be noted even the more accurately test decomposed regime has larger errors than a test decomposition in a previous study (Sellerer et al 2019), which also used a projection domain decomposition, although with higher kVp (more transmission) and larger source to image distance (less scatter). Transmission in the setup used in this work could be increased by using a higher kVp and/or adding more filtration. A new tube should be added for this and care should be taken Attenuation is given on the left y-axis and the ratio of the two is on the right y-axis in (b). The signal ratio of the two ROIs in (b) decreases as the energy increases, which is an indication that ROI 1 contains materials with higher atomic number elements than ROI 2. The contrast of the monoenergetic image (a) is defined by the energy used, with higher energies giving more weight to the soft tissue component. so that neither the dose nor the exposure time increases.
The impact of scatter was not quantified in this study, but it is a possible cause of the bias, especially if there are thick metals in the beam as discussed above. In the case of cone beam computed tomography of a head, the scatter signal can exceed the primary signal reaching the detector in the head area (Sisniega et al 2013) which would lead to severe bias in equation (5). However, panoramic imaging uses a very narrow collimation in the horizontal direction which means that the scatter is not as high as in cone beam computed tomography. Nevertheless, a strongly attenuating metallic object, such as a dental bridge (c.f. Figure 6 (a), white object on the left side of the image, on the molars), blocks virtually all direct radiation and only the overlapping scatter background remains, introducing a strong bias to the decomposition results in that area. In the future, a kernel based scatter correction method (Wang et al 2018) could be experimented with to mitigate scatter without needing a 3D reconstruction of the imaged object.
Despite the measures taken to minimize pileup (cf 2.3.3), pileup still might cause some error to the final results as charge-sharing correction increases the deadtime by an order of magnitude (Ullberg et al 2013).
Despite the aforementioned limitations, the algorithm is clearly able to separate the plastic and the skull from each other, as seen in the plot of ratios in figure 7 and in figure 6. This implies, that the limitations do not degrade neither spectral nor anatomical separation efficiency too much. A clinical study would be required to confirm if the soft and hard tissue anatomies indeed are be separated for actual human patients as well.
One implication of figure 6 is that bone mineral density measurement using a hard tissue (e.g. dentin) basis material image might be possible. It was attempted on (Langlais et al 2015) directly from the count data, but was too unstable, one possible explanation being the soft tissue background. Our algorithm takes into account the spectral response and the material decomposition yields quantitative line-integral values which might make it more suitable for this task, but further investigation is required.

Conclusion
We have built a prototype panoramic x-ray photon counting system and developed a new algorithm which enables material decomposition for panoramic imaging. By doing the calibration from the virtual panoramic images, the calibration result is enabled to be sensible despite the very low photon statistics. The decomposition of a spectral panoramic image into soft and hard tissues might be advantageous in clinical cases where the overlap of soft and hard tissues is a problem or where having the bone mineral density information might be advantageous. This should be verified in clinical trials.