A Novel and Reliable Pixel Response Correction Method (DAC-Shifting) for Spectral Photon Counting CT Imaging

1. Introduction

In the past decade, computed tomography (CT) research has been focused on replacing energy integrating detectors (EID) with energy discriminating photon counting detectors (PCD). With PCD-CT imaging, you can reach a higher spatial resolution, signal-to-noise ratio, and contrast-to-noise ratio compared to EID-CT ([1]). Furthermore, the energy discriminating ability allows you to acquire several energy bins during a scan, which allows for material decomposition, for example fat, water, and lipid concentrations ([2]), or to distinguish multiple contrast agents ([3]).

PCDs typically consist of multiple chips which each have a two-dimensional readout grid (application specific integrated circuit, ASIC) coupled to a semiconductor such as cadmium telluride (CdTe) or silicon (Si) which generates an electrical charge. The energy discrimination is implemented by thresholding the electrical charges using two or more counters. For normal operation, the first counter is set above the electrical noise of the system. Subsequent counters can be set to energy specific values, making sure that the electrical charge created by a photon is counted only when it exceeds the threshold value. One example of a PCD chip is the Medipix3RX, which consists of 256 × 256 pixels with a pixel pitch of 55 $\mathsf{\mu }$m and two threshold counters per pixel ([4]). In spectroscopic configuration, 2 × 2 pixels are combined, resulting in a 110 $\mathsf{\mu }$m pixel pitch and eight counters per pixel ([5]). The chips can be tiled together to create a large field of view (FOV). Clinical devices using this technology are the MARS Bioimaging system ([6]) and the spiral breast CT imaging system ([7])

The thresholds for the counters are set per chip and are defined as integer numbers (typical range 0 — 255 or 511), which are also known as digital-to-analogue-converter (DAC) units. The relation between DAC and actual keV threshold is close to linear and needs to be calibrated for each chip. Several methods for calibration exist, of which the following three are most common, 1) kVp scanning with a polychromatic x-ray source, where you change the DAC from high to low, searching for the setting where a pixel starts to count photons, which corresponds to the kVp ([8]); 2) x-ray fluorescence (XRF) of pure metal foils, where you search for the DAC corresponding to the characteristic K-shell electron transition energies of the metals ([8]); 3) use of monochromatic/radioactive x-ray sources, where the highest DAC with substantial counts corresponds to the source energy. With two or ideally more points within the energy range of interest, the linear relation can be determined.

While DAC settings are set per chip, the actual response varies per pixel. This means that on average, the image acquired with a certain threshold corresponds to a specific keV value, but the image will be noisy, because individual pixels represent a specific range of keV threshold values. These pixel-response variations result in ring artefacts in the reconstructed images, and therefore, projection images need to be corrected first. Basic flat field correction (FFC) is not sufficient, because the count response of the pixels is energy and flux dependent, and because of beam hardening that is present in the projection images.

One well known method to correct for the individual pixel response is the so-called signal-to-equivalent thickness calibration method (STC) ([9,10,11]). In STC, a set of homogeneous flat absorbers (for example sheets of aluminium or PMMA) are imaged with a PCD, with a fixed set of thresholds and x-ray tube settings. The non-linear relationship between absorber thickness and photon counts is subsequently fitted per pixel. To apply the STC method on imaging data of samples, you replace the counts in the projection images with the corresponding thickness. This method has the advantage that it corrects both for pixel response, and for beam-hardening. As the STC method is based on calibration data, it is highly dependent on stability of the imaging system, including the x-ray source and the detector response.

Alternative methods focus on removal of the ring artefacts either by correcting the sinogram data or the reconstructed image ([12]). Especially corrections on the sinogram data have shown to do a respectable job on removing the ring artefacts, but challenges can be present at high contrast edges, and with small objects or particles within the scanned object.

The above-described methods are mostly focused on correcting data, which is non-spectral, like with EID scans, or on correcting individual spectral threshold/bin data from PCDs. The corrections methods do not explicitly incorporate individual pixel spectral response, which is the cause of the problem in the first place. When PCDs are used to acquire multi-spectral data, with more than two thresholds, there is a risk that individual threshold correction can lead to corruption of the spectral integrity of the data.

In this work, we introduce an image post-processing method that does not require predefined calibration images. Instead, we used flat field images, acquired just before or after scan acquisition of a sample, to determine the individual pixel response in terms of DAC-to-keV at the time of scanning. In the next step, the spectral projection data was fitted per pixel and per projection using cubic-B-splines, and the counts for each pixel were corrected per keV threshold using the individual pixel response. With this method we assumed that the x-ray source and detector response was stable within the scan acquisition, while we mitigate fluctuations over longer time periods. The method allows for freedom in choosing the x-ray source settings, and the number and level of thresholds.

The purpose of this study was to describe the method, illustrate its performance, and compare the resulting reconstruction image quality with standard flat field correction, and PMMA based STC. For this purpose, we used a custom phantom with different plastics, water, and oil, and performed K-edge imaging of a phantom containing different concentrations of iodine.

2. Materials and Methods

2.1. Experimental Setup

A schematic view of the experimental setup is shown in Figure 1. A radiation shielded X-ray cabinet (EZ-Access, Metrix NDT Ltd, Loughborough, U.K.) was used, equipped with a static microfocus X-ray source with a focal spot of 5—7 $\mathsf{\mu }$m (L9421-02, Hamamatsu Photonics K.K, Hamamatsu, Japan), a sample stage with four degrees of freedom (X, Y, Z, R), and a detector stage with three degrees of freedom (X, Y, Z). A 500 $\mathsf{\mu }$m thick aluminium filter is used for all measurements to reduce the proportion of low energy photons.

For imaging, we used two Medipix3 detectors from two different manufacturers. The first was an air-cooled 2 × 2 grid layout detector with 500 $\mathsf{\mu }$m silicon (Si) sensor material biased to +150 V (Amsterdam Scientific Instruments B. V., Amsterdam, the Netherlands). The second detector had a 2 × 5 layout, water-cooling, also a 500 $\mathsf{\mu }$m Si sensor material, and the bias voltage was +150 V (WidePIX, Advacam S.R.O., Prague, Czech Republic). Both detectors have a pixel pitch of 55 $\mathsf{\mu }$m , with chips of 256 × 256, and two energy thresholds. Each chip of each detector was first calibrated using XRF measurements using pure metal foils of Cu (8.0 keV), Mo (17.4 keV), Rh (20.1 keV), and Sn (25.1 keV) ([8]);

Imaging of samples was done using step-and-shoot, meaning that the sample stage was static during acquisition of projection data. To acquire spectral data with only two thresholds per chip, we acquired multiple projection images at each rotation angle of the sample, where the thresholds in keV were changed for each projection. Before and after acquisition of projection data a set of ten open images with the same threshold settings was acquired.

2.2. Image Post-Processing (DAC-Shifting)

Our proposed image post-processing method, called DAC-shifting, consists of two major steps, see Figure 2. In step 1, the individual pixel response was modelled using the flat field data. Here we assumed an ideal setup, where the x-ray spectrum and flux are homogeneous over the entire detector. We also assumed that the XRF chip calibration assures that image acquisition with a certain threshold in keV results in comparable photon counts between the chips. With these assumptions, the ideal flat field image is a homogeneous image with a certain amount of photon counts for each threshold. For each specific threshold, the median number of photon-counts over the entire image was taken as ideal target counts per threshold. Subsequently, the individual pixel response was modelled using cubic b-spline fits through the different flat field thresholds, using the chip DAC units on the x-axis, and the pixel specific counts on the y-axis. The individual pixel fits were then used to find which DAC setting per threshold corresponds to the ideal target count. The prime outcome of this step was a look-up-table (LUT) that describes for every individual pixel which DAC unit corresponds to which keV threshold.

In step 2, the projection data was processed and corrected. For each projection angle the individual pixel data was fitted using a cubic b-spline, again using the chip DAC units of the different thresholds on the x-axis, and the pixel specific counts on the y-axis. With the LUT from step 1, the counts in each pixel and threshold were then corrected such that every pixel reflects the same keV threshold. As a result, all the projection data was corrected using the individual spectral response of the pixels. For FFC in the image reconstruction part, the ideal median counts per threshold were used instead of the actual flat fields.

It is important to note that the DAC-shifting method was performed separately for each counter. This was needed because the individual pixel response was also counter specific.

2.3. Phantom Imaging

To match the size of the detector, a small custom phantom was produced. It is a 20.0 mm diameter cylinder with a lid, both made of PMMA-c. It has seven 4.0 mm diameter holes along the length of the cylinder which hold other materials. These materials were water, sunflower-oil, PC, PEEK, POM-c, and PTFE, and centrally air, and were chosen due to their relatively similar X-ray attenuation coefficients and densities.

The phantom was imaged with both detectors, with a distance-source-detector (DSD) of 171 mm, and distance-source-object (DSO) of 123 mm. In total, 720 projections angles (0.5 degree spacing) were acquired, with for each angle five acquisitions using both threshold counters, resulting in ten threshold projections per angle in total. Exposure time per threshold was optimised to achieve similar counts for the different thresholds (Table 1).

The x-ray source was set to 50 kVp and 160 $\mathsf{\mu }$A. In each imaging session, the phantom was scanned twice. First using a static detector position, followed by an acquisition with random detector shifts between -5 and +5 pixels in the x- and y-direction, similar to ([13]). The random detector shifts were used as hardware solution to mitigate ring artefacts, which comes at the cost of more noise in the reconstructed image. For comparison between both PCDs, only data of four of the ten chips in the WidePIX detector were used. For one detector, these measurements were performed three times in a time span of 25 days. This was done to showcase the stability of the system. An overview of the measurements taken can be found in Table 2.

For image analysis, materials were segmented out using combinations of circular masks defined by centres and radii, see Figure 3. A margin (black border) was used around each selection because the materials were not perfect cylinders due to the manufacturing method and the cylinders were not perfectly aligned in the Z axis.

For the sake of illustration, Figure 4 shows an example of a CT reconstruction in all three spatial planes. The example reconstructed slice shown throughout this work was fixed at the same XY plane. This was chosen as it is a good representation of an average slice. It also contains two distinct relatively high attenuation dots of a few pixels in size in the lower right quadrant, on the edge of the PC (polycarbonate). These can be used to visually assess the reconstructed image quality. The viewing window is fixed at -1000 — 1800 HU for all slice plots, unless otherwise stated.

2.4. Image Reconstruction

For image reconstruction we used the TIGRE toolbox ([14]). Scans were reconstructed per threshold, meaning that each reconstruction represented the counts above a certain energy threshold, not an energy bin. For reconstruction we used the standard FDK (Feldkamp, Davis and Kress) option with a cosine filter. We also employed iterative FISTA (fast iterative shrinkage thresholding algorithm) reconstruction, using the FDK as prior, 100 reconstruction iterations, a hyper parameter of $2\times {10}^{-4}$, and 75 denoise iterations (‘tviter’) in every reconstruction iteration. For other parameters, the defaults were used. The chosen parameters were empirically optimised to balance the trade-off between image-sharpness and noise. All reconstructions were ${512}^3$, with isotropic voxels of 42 $\mathsf{\mu }$m.

We applied five different post-processing methods on the projection data before reconstruction. The first was standard FFC. The second method was detector based STC conversion of the counts to equivalent thicknesses of PMMA, STC-D. For STC calibration we acquired five projection images per threshold using the following thicknesses of PMMA: 0, 0.8, 1.6, 3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0, 27.0, and 30.0 mm. The median number of photon-counts over the five projections for each thickness was divided by the median counts measured at thickness 0. We subsequently modelled the individual pixel response using the global hyperbolic interpolation method from [15]. For STC-D, one STC calibration curve was fitted for every threshold, using the medians of the detector. This reconstruction step corrected for beam-hardening, but not pixel response variation. The third method was individual pixel STC conversions, where calibration curves were fitted for each pixel separately, STC-P. This reconstruction step was aimed at correcting both for beam hardening and pixel response. As fourth, the DAC-shifting method was applied to correct for the pixel response only (DAC-shifting). In a last step, the STC-D and DAC-shifting were combined to correct for beam hardening and pixel response (DAC and STC). With five post-processing options, and two reconstruction methods, each projection set was reconstructed ten times.

The imaging data of the phantom was converted to Hounsfield units (HU). HU calibration for EIDs is relatively straightforward, but this is not the case for spectral data. For PCDs with configurable energy thresholds, the HU can change with respect to the energy threshold. We chose to use the median water and air values for each energy threshold as calibration reference points. The calibration becomes:

$$HU\left(E\right)=1000\times \frac{\mu \left(E\right)-{\mu }_{water}\left(E\right)}{{\mu }_{water}\left(E\right)-{\mu }_{air}\left(E\right)}$$

(1)

Where $\mathsf{\mu }\left(E\right)$ are linear attenuation coefficients and E is the energy at which the threshold of the detector is set to. With this approach, the median water and air HU values for each threshold are fixed to 0 and -1000, respectively. This was done for every reconstruction independently.

2.5. Statistics

In this study, we compared the effect of five different options to improve image quality: reconstruction type, DAC-shifting, STC-D, STC-P, and detector motion. For presentation of the image quality, and spectral differences between the materials, the median HU over threshold was plotted, together with the 95^th percentile range of HU values present in the respective materials. For this step, cylindrical regions of interest were used, see Figure 3. The first image quality metric was the range of values for a single material within a reconstruction. We arbitrarily chose the 95^th percentile range in HU and calculated the absolute difference between the upper and lower value as the percentile size. The percentile size quantifies the precision of the system; the closer to zero it is, the better. To represent the spectral information for each material, we took the median over the different thresholds for each material within a reconstruction. To represent the overall image quality for a reconstruction, the mean percentile size over the different materials was used.

To rank the impact of each option, we evaluated their effectiveness in improving the image quality. To do so, we calculated and compared the mean relative change in percentile sizes of the reconstructions along with the interquartile ranges between each pair, without and with each parameter. These data will be named ‘effectiveness score’ in this paper and will be shown in units of percentage. For example, all FDK reconstructions versus all FISTA reconstructions. Negative values show a reduction (improvement) in percentile sizes.

In addition to this, we calculated p-values using a t-test on two related samples (without and with each parameter) to quantify the similarity of the pairs to each other. ‘This is a test for the null hypothesis that two related or repeated samples have identical average (expected) values’ [16]. A p-value below 0.01 was deemed significant.

2.6. Spectral K-Edge Experiment with an Iodine Solution

This part aimed to illustrate the capability of DAC-shifting method to enhance the spectral profile (the attenuation coefficient (${\mathrm{cm}}^{-1}$) as a function of energy (keV)) of the iodine K-edge (33.2 keV) in both 2D and 3D spectral CT imaging.

The first experiment utilised 2D spectral imaging to examine a 3D printed cuboid vessel made of polylactic acid (PLA) with a uniform cavity thickness of 6.5 mm, containing the iodine solution at a concentration of 270 mg Iodine/ml in water (Visipaque ™, Iodixanol, GE Healthcare, Finland). This setup allowed for consistent path length measurements of the iodine solution for each pixel of the detector, employing $2\times 2$ chips. Imaging was performed using the X-ray source set at 90 kVp and 89 $\mathsf{\mu }$A, with energy thresholds ranging from 8.0 to 56.0 keV in 0.5 keV increments, each with a total exposure duration of 48 seconds. A flat field correction was applied by scanning the same cuboid vessel in an empty state under identical parameters, enabling isolated assessment of iodine solution attenuation. Additionally, flat field images were used to construct a LUT essential for the initial DAC-shifting step. The outcome of this experiment focused on the effect of DAC-shifting on the spectral profile at the entire detector matrix, each discrete chip, and individual pixels.

In the subsequent experiment, a 3D-printed (polylactic acid) cylindrical phantom with a 20.0 mm diameter was utilised to perform 3D spectral CT imaging. This phantom, with nine cylindrical cavities each 4.0 mm in diameter, was filled with varying iodine concentrations (135.0, 67.50, 33.75, 16.88, and 8.45 mg/ml), water, and three empty cavities. Guided by the spectral profile of the iodine solution established in the 2D spectral imaging experiment, thresholds were selected for 3D spectral CT acquisition as presented in Table 3, specifically targeting the iodine K-edge.

The phantom was imaged with the following system configuration with detector motion and using the step-and-shoot acquisition mode:

• Tube voltage = 90 kVp
• Tube current = 89 $\mathsf{\mu }$A
• DSO (distance source object distance) = 123.0 mm
• DSD (distance source detector) = 163.0 mm
• Projections = 720
• Rotation angle = ${360}^o$
• Reconstructed voxel size = 0.0415 ${\mathrm{mm}}^3$

Ten flat field images were acquired before and after the projection images. CT reconstructions were then performed as described in Section 2.4 using both reconstruction algorithms, without and with DAC-shifting. Cylindrical regions of interest (volume = 79 ${\mathrm{mm}}^3$) as shown in Figure 5 were used to evaluate each iodine concentration.

The quantification of different iodine concentrations was achieved by examining the Hounsfield units at 33.0 keV, the energy level closest to the K-edge peak. A linear regression model ($y=mx$) was used to evaluate the correlation between iodine concentrations and CT numbers. The y-intercept was set to the HU of water (zero). Further, the ${\mathrm{R}}^2$ values of without and with DAC-shifting methods at 33.0 keV were compared using Fisher’s r-to-z transformation, calculating the z-score and corresponding p-value to determine the significance of the difference between the two correlation coefficients.

3. Results

3.1. Out-of-the-Box Imaging

To set the reference, we start with a CT reconstruction showing the out-of-the-box (OOTB) experience (upper row of Figure 6). Apart from manufacturer specific procedures, the only calibration done at this stage was a basic energy calibration per chip using XRF. This only ensured that the chips were aligned in terms of energy to first order.

The spectral profile of the different materials can be seen in Figure 7. Note that despite the increasing exposure time used for increasing thresholds, there was a substantial increase in percentile size with increasing thresholds. Also note that the 95^th percentile bands overlap each other for most materials.

3.2. Effectiveness Scores

In Table 4 an overview of the effectiveness scores for the investigated options on the different datasets can be seen. In the following sections each option is separately illustrated.

3.3. STC-D Beam Hardening Correction and Detector Motion

The effect of the STC-D beam hardening correction is visualised in Figure 6. From the profile lines in the figure, especially at the 5 keV threshold, clearly there was beam hardening present in OOTB, which was corrected using STC-D. The overall effect of STC-D on the effectiveness score was limited (Table 4), as the beam hardening was present mainly in the PMMA bulk and the PTFE, but not in the other materials.

The detector motion method effectively translated the ring artefacts caused by pixel response differences into random appearing noise in the reconstructions (Figure 8). As the method changes the appearance of the image, but not the value range within the materials, the mean percentile sizes were not substantially affected most of the time (Figure 9). The combined effectiveness score for detector motion limited to -8.7 % (IQR -13.0 — -1.5 %), $P<0.001$.

3.4. STC-P Beam Hardening and Pixel Response Correction

We observed substantial differences in the effectiveness of STC-P with the different datasets (Figure 10). The effect of STC-P on two different datasets is visualised in Figure 11. The corresponding effectiveness scores varied in range from -10.0 — -43.7 % for the ASI data, and even a worse image quality for the ADVACAM-1 data (effectiveness score +45.5 % (+39.2 : +52.3), Table 4). The two datasets that were acquired close to the date of acquiring the STC-D calibration measurements (ASI-2 and ASI-3) seemed to benefit more from STC-D than the other datasets.

3.5. DAC-Shifting

The DAC-shifting method was focused on correcting the individual pixel response, using a LUT based on the flat field images of the day. An example of a LUT for a specific threshold and its variation over different datasets can be seen in Figure 12. In this example, the 95 % range of DAC settings for individual pixels that correspond to 14 keV stretches from 73 — 105 DAC units. The stability over different datasets seems random and the root-mean-square (RMS) standard deviation of Figure 12 c corresponds to 2.33 DAC units, which typically corresponds to 0.47 keV.

In Figure 13, we illustrated the effect of DAC-shifting in projections and sinograms at two thresholds and compared it to STC-P. Note the significant reduction in vertical stripes in the sinogram images, especially at the higher threshold. In addition, individual chips can be distinguished in the OOTB and STC-P projections, but not in the DAC-shifting projections.

With DAC-shifting, ring artefacts were substantially reduced, without degradation of spatial resolution or introduction of obvious artefacts (Figure 14). DAC-shifting also resulted in comparable image quality between different datasets (Figure 15). The absolute mean percentile sizes with DAC-shifting are comparable between the datasets, even when the absolute numbers without DAC-shifting vary (for example, ASI-1 vs. ASI-3). As a result, the effectiveness scores were substantial (mean -47.4 %) and varying based on the initial image quality (-40.7 — -55.6 %), $p<0.001$.

3.6. Reconstruction Methods

So far, we have visualised reconstructions using only the FDK method. The comparison between the iterative reconstruction algorithm, FISTA, and FDK is shown in Figure 16. The FISTA method as used in this work was very effective in reducing the noise in the image, with minimal compromise for the image sharpness. The noise reduction clearly translated into reduced mean percentile differences. As a result, the effectiveness score was -61.6 % (IQR -66.0 — -59.1 %), $P<0.001$.

3.7. Optimal Settings

The best combination of options was detector motion, STC-D, DAC-shifting, and FISTA reconstruction. Even though the effectiveness score for detector motion was close to zero, it complements the FISTA reconstruction. This is because detector motion smears the rings into noise, which can effectively be reduced by the FISTA reconstruction. An illustration of the image improvement compared to the out-of-the-box setting using FISTA can be seen in Figure 17. The spectral data of the best parameter combination can be seen in Figure 18.

3.8. Evaluation of DAC Shifting - Iodine K-Edge Imaging

The capability of the DAC-shifting method to enhance the sharpness of the iodine K-edge using 2D spectral imaging is illustrated in Figure 19. The spectral profiles show that attenuation increases as the energy approaches the K-edge and decreases after passing it. Without DAC-shifting, individual pixel response variations led to a relatively wide 95 % data range. In contrast, the DAC-shifting method improved homogeneity between individual chips (see attenuation images) and pixels (see 95 % data range), resulting in a sharper attenuation peak near the K-edge.

In 3D spectral iodine imaging, the spectral profiles, depicted in CT number (HU) as a function of energy, still show increased CT numbers around the K-edge, but less pronounced than in the 2D imaging (Figure 20). Note that the highest concentration in the 3D imaging was half of the concentration in the 2D imaging. The effect of DAC-shifting was again mainly visible in a reduction of the 95 % data range, especially at the higher thresholds. This resulted in less overlap between the 95 % data ranges in the lower iodine concentrations, as can be appreciated from the linear relationships depicted in Figure 20 Comparison of ${\mathrm{R}}^2$ values demonstrating a highly significant difference (p <0.001) between without DAC-shifting and with DAC-shifting methods.

4. Discussion

In this paper, we presented a novel post-processing method to deal with individual pixel response in spectral imaging with PCDs. We have shown how the DAC-shifting method has a positive impact on the image quality, also in relation to existing methods. We have also shown that DAC-shifting is robust against long term variations in the imaging system, while STC-P did not provide consistent image-quality improvements over time.

When comparing the DAC-shifting method with the STC-P method, we can see that at best, with the most optimal settings, they perform almost comparably. In the ASI datasets with motion, the mean percentile sizes for STC-P with FISTA were 68, 57 and 53 HU for ASI-1, ASI-2 and ASI-3, respectively (Figure 10). For DAC-shifting in combination with STC-D and FISTA reconstruction, the mean percentile sizes were 50, 50 and 49 HU (Figure 15). In general, the results on DAC-shifting were very stable when comparing the effect of the different post-processing methods on the different datasets, even when comparing the two detectors (Figure 15). For STC-P, the changes in mean percentile size were much less stable, and it appears that it only worked for ASI-2 and ASI-3, which were acquired within four days of the STC calibration measurements. These findings exemplify the instability of this imaging system, this is observed in previous work using a Medipix3 based system [17]. This instability was further illustrated and quantified in the LUT variation example in Figure 12 where the RMS of the standard deviation map was 2.33 DAC units, which corresponds to a 0.47 keV standard deviation in this system. In other words, within the data acquired in a four week window, individual pixel responses at 14.0 keV can vary between -1.0 — +1.0 keV. Given the random distribution of the variation in the example image, the variation can be attributed to the detector, not the x-ray source.

Even though several papers are available on STC-P, none of them describe the relation between system stability and effectiveness ([9,10,11,15]). This makes a direct comparison to other papers challenging. The most comparable work in terms of imaging with Medipix3 and low threshold energies is the work by Ronaldsen et al. using the MARS-CT scanner in 2012 [2]. The spectral profile on sunflower oil obtained in our work (Figure 18) seems directly comparable to their work, with slowly increasing CT numbers above 15 keV threshold. Unfortunately, they only mention a standard deviation uncertainty of ~10 HU, without explicitly showing the data. They make use of a combined wavelet-Fourier filtering method for ring artefact correction [18], but the effect of this is not separately illustrated. They also imaged iodine, at a concentration of 0.01 mol/l, lower than our lowest concentration (0.033 mol/l). However, they only measured thresholds of 30.9 and 36.2 keV around the K-edge, for which the K-edge was not directly visualised.

He et al. performed more comparable K-edge imaging of iodine, both in 2D and 3D [6]. They showed similar graphs regarding the K-edge in terms of attenuation over energy (Figure 19). In their and our work, the ability to precisely resolve the K-edge is limited by the effective spectral resolution of the detector. The spectral profile appears to be between the theoretical mass attenuation example with a 1.0 and 3.0 keV ${\mathsf{\sigma }}_E$ (energy resolutions) provided by Ge et al. [6]. In their example, they show that a limited energy resolution convoluted with the sharp K-edge, which shifts the attenuation peak to a slightly higher energy.

When calibrating the Medipix3 detector, an equalisation procedure minimises the variation in pixel response at a specific threshold. In our setup, we focus our research on low keV thresholds. The Advacam detector was equalised at approximately 8.0 keV, whereas the ASI detector was configured using the noise edge, which is approximately 4.0 keV. Due to threshold dispersion, where the pixel responses diverge as the threshold is increased, the DAC-shifting improvement was more significant at the higher thresholds, achieving comparable image quality between 5.0 and 20.0 keV (Figure 17). Alternative methods, like loading specific equalisation maps for each threshold used could help increase the image quality of the raw data. Currently, this step would increase the acquisition time substantially, as loading of the equalisation maps takes time.

Due to the instability of our imaging system, the STC-P method did not provide consistent image quality improvement. To make use of the advantage of STC in correcting for beam-hardening, we introduced the STC-D method. To the best of our knowledge, this approach was not used before. In our setting, the advantage was that STC-D can be combined with the DAC-shifting method. The method is simple, reliable and fast, because it can be applied to the entire stack of projection data at once and is less sensitive to variations over time. Of course, the effectiveness of STC-D beam hardening correction is dependent on the materials used for calibration, PMMA in our case, and the sample that is imaged. For the phantom with plastic inserts it worked well, but for the iodine acquisitions, it showed less effective (data not shown), this was at least partly due to the minimal beam hardening at the higher energies.

In recent years, there has been focus on improving spectral image quality by utilising dedicated spectral reconstruction methods with material decomposition [19,20,21], and deep-learning based reconstruction for noise reduction [22,23], or combinations thereof [24]. While these novel methods all seem very promising, they could still use the advantage of the presented DAC-shifting method. The method corrects the projection data, keeping the spectral integrity, without changing the units, and homogenises data between the different chips and pixels. In this way, the input data for reconstruction is improved, providing improvements for downstream steps like material decomposition. The challenge of our method is that it adds measurement and processing time. It took between 10 and 12 minutes to apply DAC-shifting on a modern 64 core computer (AMD Ryzen Threadripper 3990X) utilising Python and parallel processing on 60 of the cores for most of the steps. To speed this up, there are many options. One of which is a neural network for the DAC-shifting method to improve processing speed.

The presented DAC-shifting method has some limitations. First of all, the method is based on performing cubic b-spline fits through the spectral data, which means that it can only be used on spectral data, and not for individual single-threshold data. The cubic b-spline approach does not utilise any knowledge about the underlying physics. Ideally we would fit a function which is related to the underlying physics with a data based uncertainty estimation to constrain the fits. Secondly, at the edges of the spectral data, the lowest and highest threshold, approximately half of the pixels were corrected using extrapolation of the data, which is less accurate than interpolation due to the cubic b-spline approach. This effect was visible in Figure 18, where the 95 % data range increased at the highest threshold plotted. Note that in our work we already used one extra threshold acquisition (Table 1 $\#5$, threshold 0 = 25.0 keV; threshold 1 = 30.0 keV), to mitigate this issue. This data was only used to guide the spline fitting, but the DAC-shifted reconstructions with these thresholds did not result in improved image quality. This is due to the extrapolation (data not shown). Thirdly, since each counter corresponds to different transistors, and therefore each counter has their individual response patterns, the DAC-shifting method was applied per counter. This means that the approach is currently not applicable to PCDs that have many more counters per pixel where the entire spectral profile can be acquired at once. Further research is needed to address application in these detectors, a first approach was taken in [17]. Finally, we reconstructed our data as spectral thresholds, not spectral bands (energy bins), which are often presented in other papers. However, the DAC-shifting method is applicable in both cases. When spectral bands are needed, the DAC-shifting method can be applied first on the threshold data, and subsequently spectral bands can be calculated by subtracting two thresholds.

5. Conclusions

Within the context of improving reconstruction quality of our system, we have presented a novel post-processing method that minimises individual pixel response in spectral photon counting cone-beam CT. The method was shown to be consistently effective in correcting the sinogram data, which resulted in reduced ring artefacts and improved image quality. This was supported with the iodine concentration measurements. It does not appear to introduce artefacts and maintains the scale and units of the data. DAC-shifting was also shown to provide consistent image quality, even with variable quality input data.