Ring artifact suppression in X-ray computed tomography using a simple, pixel-wise response correction

We present a pixel-specific, measurement-driven correction that effectively reduces errors in detector response that give rise to the ring artifacts commonly seen in X-ray computed tomography (CT) scans. This correction is easy to implement, suppresses CT artifacts significantly, and is effective enough for use with both absorption and phase contrast imaging. It can be used as a standalone correction or in conjunction with existing ring artifact removal algorithms to further improve image quality. We validate this method using two X-ray CT data sets acquired using monochromatic sources, showing post-correction signal-to-noise increases of up to 55%, and we define an image quality metric to use specifically for the assessment of ring artifact suppression. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

used aluminum filters to attenuate the beam and took the assumption that the spatial variations in response are primarily caused by changes in thickness across the scintillator. The algorithm presented here makes no assumptions regarding the cause of the detector's spatial variations and uses a method that is sample-independent; it has been tested here for the specific case of a monochromatic source. Here, we focus on the problems associated with detector defects and changes in the optical system, presenting a simple, pixel-wise detector calibration using hundreds of data points for each position on the detector, rather than the standard two-point flat-field calibration. This method typically requires only a single image sequence for each experimental setup and does not require any information about the physical cause of the detector's response variations. We assume, for the sake of simplicity, that any correlations between adjacent pixels can be considered negligible for this purpose, and as our results will show in section 7, this is found to be a reasonable assumption. Taking the detector point spread function into consideration in our reconstruction would likely further improve the results, but doing so would require a considerably more sophisticated approach that is beyond the scope of our study.
Our approach can be used for conventional X-ray CT, and indeed would even be useful just for a more accurate representation of pixel values in 2D imaging, but it was primarily motivated for the case of phase-contrast X-ray CT [22]. Phase contrast imaging is a relatively recent development, which converts the phase shifts imparted to the X-ray wavefield by the sample into intensity variations that can reveal soft tissue features. The simplest way to achieve this additional soft-tissue contrast is through the introduction of a distance (e.g. > 1 m) between the sample and detector, illuminating the sample with sufficiently coherent X-ray radiation [23,24]. Since phase-contrast X-ray imaging can capture soft tissue structures, we have recently applied the method to image the brain in situ, finding that the phase contrast associated with grey and white matter is orders of magnitude less than the attenuation contrast associated with the skull [25]. The correction presented here is therefore of particular interest in phase contrast imaging, which can reveal weakly absorbing objects whose contrast is low enough to be comparable to the contrast of ring artifacts. This requires correction at a level of contrast that is typically buried within the noise of conventional CT. For this reason, post hoc methods that generally perform quite well may not be as effective for low-noise CT data, such as can be achieved with phase retrieval. A comparison between some common algorithms and the one presented here can be found in the Appendix.
Phase retrieval removes the effects of phase contrast to provide quantitative reconstruction of the attenuation coefficients, while massively suppressing high frequency noise in the reconstruction [25][26][27][28]. This can enable the radiation dose to be greatly reduced [28] and enables brain tissue to be resolved in situ in the heads of small animals [25], as seen below. While the noise suppression of phase retrieval will reduce the noise associated with ring artifacts, it can also make ring artifacts seem more prominent than for absorption contrast imaging by greatly suppressing the random noise that can hide ring artifacts [25].
During the review of this manuscript, a paper was published that takes a similar, though in some ways less general, approach to ours [29]. In that work, the authors show their pixel-wise approach to be effective for absorption contrast imaging based on the more uniform illumination of a conventional X-ray source.

Primary, residual, and external factor corrections
The phenomena that alter the detector response occur over a range of different time scales. We consider here effects over three distinct time frames: 1) Long-term effects that have accumulated over the lifetime of the optical system components. These describe the initial state of the system at the start of an experiment and persist throughout. An optical system will almost always have some number of pre-existing defects at the beginning of an experiment. For example, the scintillator may have scratches on it or other thickness variations. There may be oil or mineral residue on the surface of the optics from handling and/or moisture deposits, and dust can collect on any part of the system, particularly the scintillator, which is often exposed to the environment. These imperfections are the easiest to correct and can be characterized at the beginning of an experiment, since they will persist within the system throughout the duration of the experiment. 2) Experiment-scale effects that arise over the course of the time during which an experiment takes place (e.g. dust accumulation on the scintillator [30] or variations in dark current with temperature changes). Effects on this scale can be accounted for separately. 3) Short-term effects that change during a single CT acquisition. This is most commonly seen with beam injections in synchrotron data that occur on a timescale much smaller than the acquisition time but much longer than the exposure time of the individual projections. In addition, dust and other imperfections can also accumulate over the course of a CT acquisition. The latter is not usually a problem for most CTs lasting on the order of a few minutes to an hour or two; however, it may be a problem for longer CT acquisitions on the scale of hours to days.
Furthermore, it is important to differentiate between phenomena that are directly related to the detector system (e.g. settled dust or scratches) and those that are external to the detector system (e.g. fluctuations in the X-ray source or polychromaticity). The primary purpose of this work is to correct for those effects that arise within the detector system itself, in a very simple way. Every experiment will also have different external conditions that must also be considered. Here, the intensity variations caused by the synchrotron beam injections are non-negligible and must also be taken into account, so we have included a correction for them along with the corrections that are specific to the detector system. As we have developed the primary correction here using a monochromatic source, and the detector response can vary with energy, the matter of polychromaticity is left as a subject of future work.
The full correction described herein consists of three distinct corrections. The first, which we will refer to as the 'primary' correction and which is the main focus of this work, accounts for both pixel-to-pixel variations in the detector response as well as any pre-existing imperfections that affect the detector response (e.g. scratches and/or dust on the scintillator); this correction accounts for the state of the detector system at the start of the experiments (see (1) above). The second part is a 'residual' correction and accounts for any changes to the detector system that occur over the course of an experiment (2), which can be several days after the primary correction has been defined. Note that, while this correction cannot account for changes in response that occur over the course of a single CT, these types of short-term changes in the detector system are less common than those that occur over the longer time frame of the experiments and often do not pose an issue over the typical CT acquisition time. The third and final correction is a per-projection scaling for source-intensity variations, which occur on a timescale shorter than a typical CT (3).

Pixel-wise mapping of the initial detector response
The primary correction can be broken into two main parts -first, a spatial mapping of the detector response, followed by an application of that mapping to experimental data. In the following sections, we break these parts further into eight basic steps, labeled (1) -(8) below. These steps are summarized in the flow chart of Fig. 1. To determine the spatial variations in response across the detector that result in ring artifacts, we only need to know how the intensity measured by the optical system at each pixel varies from that which is incident upon it. We take a series of measurements acquired across a range of intensities covering the bulk of the dynamic range of the detector. It is possible to use filters to attenuate the beam [21], however spatial variations in density and thickness of the filters can add artifacts back into the images.
Step 1 -To avoid introducing new artifacts, we start by simply sweeping the detector through the X-ray beam while acquiring a sequence of images (see Fig. 2). Since the beam intensity in the absence of any sample is typically peaked at the center, rolling off outward away from   4) and (5) The measured intensity is fit as a function of the 'true' intensity for every pixel, yielding detector gain and dark-current offset maps. (6) For a CT data set, projections are corrected using these maps, after dark-current correction. (7) The mean of the flat-field images acquired with the CT data is dark-corrected, and the resultant image is smoothed using the same smoothing kernel as that used in step 3. (8) The corrected projection images are then flat-field corrected using the smoothed flat-field image, if necessary with scaling for beam injections. *The residual correction (7) and beam injection scaling (R P /R F ) may be omitted if not required.
the beam, a wide range of intensity measurements can be obtained. For a beam size that is larger than the detector, a single sweep is usually sufficient; when the detector is larger, more than one sweep with the beam offset from the center of the detector may be necessary to ensure coverage across the full range of intensities for every pixel. By eliminating the need for absorbing materials, we can measure the exact response for a given input intensity, thus avoiding any possible artifacts that may be introduced by structural imperfections in the absorber. If, however, the intensity profile of the beam is insufficient to sample the full range of intensities required, then filters can be added to extend the intensity range acquired; the structural imperfections would not introduce additional artifacts when used in this way, since they would be shifted across the extent of the detector in the sweep direction during acquisition. It should be noted that the same measurements can be acquired in a number of ways. While, for this work, we endeavored to eliminate the absorber entirely, an image sequence can also be acquired by rotating an absorber such that the attenuation varies with the thickness of absorber along the line of sight [14] while preventing density non-uniformities in a static absorber from introducing new ring artifacts. For a conventional X-ray source, the same effect can be achieved by varying the source current.  The experimental set-up, shown for an incident beam with a parabolic intensity profile. When sweeping vertically, the detector is positioned high enough above the beam to achieve the desired minimum number of incident counts. A sequence of images is then acquired while the detector is translated downward through the beam to the equivalent position below the beam. Three images are shown near the middle of the sequence, along with intensity profiles taken vertically through the center of the images. The extreme ends of the sequence (not shown) contain only traces of the edge of the incident beam. Note that the direction of transverse translation is not important; it is equally valid to acquire an image sequence while translating the detector horizontally through the beam.
Image sequences were acquired on beamline BL20B2 at the SPring-8 synchrotron in Hyogo, Japan using a 2048 × 2048 ORCA Flash 4.0 digital sCMOS camera (C11440-22C by Hamamatsu) with a 25 µm thick gadolinium oxysulfide (GOS) scintillator coupled with a tandem lens system, giving an effective pixel size of 15.1 µm. Three beam sweeps were required for full coverage across the detector. 270 images were acquired in each sequence, with an exposure time of 100 ms each. Sequences were recorded at three positions such that the beam was centered on the horizontal left, center, and right side of the detector, which was swept vertically.
Step 2 -Once the full beam sweep sequence was acquired, each image was dark-corrected using the mean image D(i, j) of the dark-current images acquired at the time of the sequence, and then the images were stacked into a volume I(i, j, k), where i and j are the spatial dimensions within the k images in the stack. Figure 3(a) shows the intensity through the combined image stack for a single pixel.
Step 3 -A separate volume was created wherein each image I(i, j) of the volume I(i, j, k) was smoothed with a Gaussian filter, yielding an estimate I s (i, j, k) of the true intensity profile of the incident X-ray beam. The smoothing kernel radius was chosen to be large enough to eliminate the pixel-to-pixel variations resulting from the non-uniform response, while being small enough to maintain the underlying shape of the beam intensity profile. We found that a kernel with a standard deviation of 50 pixels provided a suitable blurring function for our experimental conditions.
Step 4 -The 'true' beam intensity was plotted against the measured intensity to determine the order of polynomial, for smoothed counts as a function of measured counts, required to model the response. This is shown in Fig. 3(b). Since the two intensities come from the same measurement, this resultant curve is necessarily linear, with a slope close to unity. The 'true' intensity I s (i, j, k) was fit for each pixel (i, j) as a linear function of the measured intensity, Step 5 -The coefficient arrays α = α(i, j) and β = β(i, j) from this calibration (specified in bold for clarity) can then be used to determine the 'true' intensity of the beam I s (i, j) incident on the detector, given the measured intensity I(i, j) for a given projection image, regardless of the sample. It should be noted that this will correct for spatial variations in the detector response and offset, which are those primarily responsible for ring artifacts; however, since I s is estimated by smoothing the data itself, this method cannot account for any large-scale non-linearities (e.g. due to higher-order harmonics or polychromaticity); these must be accounted for separately. However, these effects are expected to make a smaller contribution to ring artifacts than the small-scale, pixel-to-pixel variations caused by imperfections in the optical system.  Fig. 3. a) The counts measured at a single pixel in a stack of images taken while the detector was swept vertically through the beam three times; the beam was centered first on the horizontal left, then center, and finally the right side of the detector. Note that the highest peak for this particular pixel corresponds to the sweep across the center of the detector, with the next highest peak corresponding to the right, indicating that this pixel is located just to the right of the horizontal center of the detector. b) The 'true' (i.e. smoothed) counts for the same pixel as in (a) as a function of the measured counts.
The coefficient maps, α and β, are shown in Figs. 4(a) and 4(b), respectively. Note that α corresponds to the detector gain, while β maps the intercepts of the fits and hence an offset to the dark current. A distinct line can be seen clearly across the upper left of Fig. 4(b), and to a lesser degree in Fig. 4(a), possibly -though not necessarily -due to a scratch on the scintillator, as per the assumption of scintillator thickness variations of Vågberg et al. (2017) [21]. Our correction does not require any assumption about the cause of these response variations, so while effects resulting from thickness variations in the scintillator may be present, any other phenomena affecting the response would be represented as well.
To test our method under different experimental conditions, a calibration was also performed for a different detector at the Imaging and Medical Beamline (IMBL) of the Australian Synchrotron. Image sequences were acquired using a 2560 × 2160 pco.edge 5.5 sCMOS camera with a tandem lens configuration and a 25 µm thick GOS scintillator, giving an effective pixel size of 16.2 µm, at a beam energy of 25 keV. Due to the relatively small beam height and narrow vertical slit aperture, five separate image sequences of 300 exposures each were acquired while sweeping the detector across the beam horizontally, rather than vertically, with each sweep centered at a different vertical position on the detector. The images in each sequence were segmented into five strips of equal height, covering the full width of the detector. The strips from each sweep corresponding to the region of peak intensity were combined to form a single final image sequence. This was done to provide a single beam sweep sequence with full coverage across the intensity range at every pixel, while removing the vertical slit edges present in the individual sequences. From this combined sequence, smoothed and unsmoothed volumes were created and processed following the same method as above, yielding the gain and dark-current offset correction maps shown in Figs. 4(c) and 4(d), respectively.

Applying the primary correction
Step 6 -To implement the correction for a CT sequence, Eq. 1 is applied using the coefficient arrays from the beam sweep and replacing the beam sweep volume I(i, j, k) with the projection images P(i, j, k) and using the dark current acquired during the CT sequence, yielding the corrected projections: (2)

Residual correction
Step 7 -A new flat-field image is created from the mean of the flat-field images acquired with the CT sequence, smoothing the result using the same smoothing parameters as those used to create the beam sweep volume I s . If the state of the optical system has not changed between the time of the beam sweep and that of the CT, then this smoothed flat-field image can be used for flat correction. However, if the state has changed, then an additional correction may be required. This is done by applying the same correction above (Eq. 2) to the unsmoothed mean flat-field image, with the residuals calculated as the difference between the corrected and uncorrected images. These residuals, which represent the changes in the system over the course of the experiment, are then added to the smoothed flat-field image used for the flat correction.
It should be noted that the residual correction is dependent on the amount of noise within the data set. There will likely be unwanted signal within the residuals (e.g. zingers, pixels where the intensity is far from the intensity in surrounding pixels), and when the noise level is low, adding these residuals to the smoothed flat-field image can introduce ring artifacts. In this case, it is important to filter the residuals to ensure that only those that are persistent throughout the CT are included. This can be done by limiting the inclusion of residuals only beyond a certain tolerance (we use those > 3σ from the mean). When the data set is noisy, this filtering is less important and can even remove residuals that should be included.
Step 8 -Finally, the corrected projection images are flat-field corrected using the flat-field image created in Step 7. It can often be useful to correct for changes in the beam intensity over the course of the CT acquisition due to beam injections by customizing the flat-field image for each projection. One common way to do this, and the method used here, is to scale the smoothed flat-field image by the ratio of a reference region in the projections outside the sample to that same region of the mean flat-field image.

Experiments
During the synchrotron experiments described in the previous section, propagation-based, phasecontrast CTs were acquired of biological samples under the same experimental conditions as those for the beam sweeps. A CT sequence was acquired at 24 keV at SPring-8 of a scavenged head from a New Zealand White rabbit kitten born at 30 days gestational age (GA; term ∼32 days), suspended in agarose. 1800 projections were acquired at a 5 m sample-to-detector propagation distance over 180 • , with an exposure time of 100 ms per projection. A second phase-contrast CT sequence was acquired at the Australian Synchrotron of the lungs of a New Zealand White rabbit kitten, also born at 30 days GA, at 25 keV using a MICROFIL contrast agent. 3600 projections were acquired at 2 m propagation over 180 • , with an exposure time of 200 ms per projection. To evaluate the effectiveness of our correction, each data set was processed twice, once with traditional dark-current and flat-field corrections and once with the pixel-wise correction detailed in this paper.
Phase retrieval was performed using the two-material algorithm derived by Beltran et al. (2010) [27] from the single-material algorithm of Paganin et al. (2002) [26] and described for CT by Croton et al. (2018) [25]. For the head data set, the bone and soft tissue interface was retrieved, while for the lung data set, phase retrieval was performed with respect to the MICROFIL/tissue interface. In total, four volumes were reconstructed for each data set -before/after correction without phase retrieval and before/after correction with phase retrieval.

Results
Several iterations of a reconstructed, phase-retrieved slice from the rabbit kitten data set are shown in Fig. 5, with different steps of the full ring artifact correction applied in each. Images in the top row were created using (a) a standard dark-current and flat-field correction, (b) applying only the beam injection correction, (c) only the primary response correction, (d) the primary response and residual corrections, and (e) with all three of the primary response, residual, and beam injection corrections applied. The bottom row shows (f) the sinogram for the standard dark current and flat-field correction and (g) -(j) the differences between (f) and the sinograms for each of the corrections applied in the top row. The intensities of the reconstructed tomograms are all scaled to the same palette, and the difference images are scaled independently but with the same size intensity interval for comparison. The ring artifacts in Fig. 5(a) manifest as vertical stripes in the sinogram in Fig. 5(f), however they are imperceptible in the latter image, since they are overwhelmed by the signal from the skull. This makes many of the sinogram-filtering techniques less effective. The difference images of Fig. 5(h) -(j) clearly show the rings that have been removed by the response correction described herein as well as the temporal variations in intensity, seen as horizontal stripes in panel (j), due to the beam injections over the~3-minute duration of the CT acquisition. In addition, because sample features appear in Fig. 5(j), it can be seen that we are not just correcting these artifacts, but we are also correcting the intensity of the images. While stripes due to ring structures can be seen clearly in the difference sinogram with only the beam injection applied (g), the contrast between these artifacts is quite low relative to the other corrections. Indeed, in the reconstructed tomogram (b), very little difference can be seen compared to the standard corrections (a). This correction appears to primarily affect the contrast of the image rather than the artifacts themselves. This is not surprising, given that the changes due to the beam injections (horizontal lines in (g)) are orthogonal to the time axis. The most striking changes with respect to the ring artifacts can be seen for the primary response correction in (c) and (h). The ring artifacts are greatly reduced; however, a bright ring artifact remains (yellow arrow), most likely as a result of changes in the time frame between detector characterization and the CT scan. This effect is not common throughout the data set, and this particular tomogram was explicitly selected to demonstrate the effects of changes to the optical system between the time that the beam sweep for the primary response correction was acquired and the time that the CT was acquired. Once the residual correction is included (d), this artifact is removed completely. Frames (e) and (j) show the result with all three of the primary response, residual, and beam injections applied.
Figures 6(a) -6(l) show reconstructed slices of the rabbit kitten head and lung CTs, both before and after phase retrieval, for standard dark-current and flat-field corrections and for the pixel-wise response correction described in section 4. The phase contrast tomograms in the top row are included to demonstrate how small the effect of the ring artifacts is without phase retrieval (similar to absorption contrast) compared to when phase retrieval is applied. The contrast of the ring artifacts is only slightly above the noise, and they become substantially more prominent once that high-frequency noise is suppressed with phase retrieval. Slices were also reconstructed for each data set, post-phase-retrieval, of a sample-free region containing only the sample container filled with agarose. In each of the data sets, the signal-to-noise ratio (SNR), measured as the ratio of the mean to the standard deviation (SNR = µ/σ) [31] within the region of interest, increased significantly in the area immediately surrounding the center of rotation (COR) after the correction was applied. We refer the reader to references [28,[32][33][34][35] for further details regarding the SNR boost associated with phase retrieval. The artifacts are strongest near the COR and fall off with increasing radius, since incorrect pixel values are spread across increasingly larger circumferences. In the innermost region, within a radius of 100 pixels from the COR, SNR increases of 38% and 39% were achieved with the correction for the phase-retrieved head and lung images shown in Fig. 6, respectively, and increases of up to 55% were seen across the full sample-free volumes.
To better quantify the improvement, we define an image quality metric, the ring artifact suppression percentage (RASP), as the percentage reduction in the standard deviation (σ) of the azimuthally averaged radial profile from the COR in the sample-free region: Here, c and u denote the corrected and uncorrected images, respectively, N is the total number of  All data are shown without (first and third columns) and with (second and fourth columns) the correction detailed in this paper. Top row: Phase contrast tomograms of rabbit kitten head and lungs, no phase retrieval. Middle row: Tomograms of the same rabbit kitten head and lungs, after two-material phase retrieval [25,27]. Bottom row: Phase-retrieved tomograms of sample-free regions from rabbit kitten head and lung data sets. First and third columns: Standard dark-current and flat-field correction. Second and fourth columns: Beam sweep response correction. radial bins, x j is the azimuthally-averaged value within each bin, and µ j is the mean value in each bin, determined by smoothing the radial profile with a median filter.

Conclusions
We have presented a simple correction to account for small spatial variations in detector response that is both easy to implement and highly effective at removing CT ring artifacts. This method is applicable for both absorption contrast and phase contrast imaging, and artifacts are suppressed to a level low enough to enhance even very low-contrast features, such as the soft-tissue boundaries a b d c Fig. 7. a) Azimuthally-averaged radial profiles for the sample-free region in Figs. 6(i) and 6(j). The decrease in the variance between the uncorrected and corrected images can be clearly seen, particularly at smaller radii, where the effects of the ring artifacts are strongest. b) The reconstruction region used for averaging (Hamamatsu ORCA Flash 4.0). Since the CT was acquired over 180 • , the artifacts for individual pixels occur only in the top or the bottom half of the image. c) Azimuthally-averaged radial profiles for the sample-free region in Figs. 6(k) and 6(l). d) The reconstruction region used for averaging (pco.edge 5.5). The images in (b) and (d) are the uncorrected tomograms.
within the brain. This correction requires just a single series of images to be acquired while sweeping the detector through the X-ray beam, and the data need only be acquired once for each experiment with a common detector configuration throughout. In addition, the algorithm should be applicable to both laboratory and synchrotron experiments alike, with some modifications required for a polychromatic source. This correction is effective on its own but can also easily be combined with other existing ring artifact removal methods for further improvement. We have also defined a ring artifact suppression metric that can be used to assess the quality of any ring artifact removal technique, and we have used this metric to quantify the improvements achieved with our method. Figure 8 shows a comparison between several different ring artifact removal methods applied to phase contrast imaging data. These methods were implemented using the TomoPy Python package [36] in conjunction with the ASTRA Toolbox [37]. Every method was implemented several times, each time with its parameters varied. Between three and fifteen corrections were created for each method, depending on the number of parameters required and how quickly each successive attempt converged toward a cleanest reconstruction; the most improved results are shown. Methods employed included a smoothing filter (size=15), a combined wavelet-Fourier filtering method [8] using a db5 filter (level=4, sigma=5), a Haar filter (level=1, sigma=4), and a sym5 filter (level=4, sigma=2), and Titarenko's method (alpha=3) [11].  Fig. 8. A phase-retrieved tomogram of a rabbit kitten brain in situ, with five different ring artifact correction methods applied, as labeled. The image on the upper left is the uncorrected image, with only a standard dark-current and flat-field correction applied, before phase retrieval. The image on the lower right has been corrected following the procedure outlined in this paper.