{'fraction_non_alphanumeric': 0.03430034129692833, 'fraction_numerical': 0.008617747440273037, 'mean_word_length': 4.656611969111969, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low-rank flat-field correction for artifact reduction in spectral computed tomography

Spectral computed tomography has received considerable interest in recent years since spectral measurements contain much richer information about the object of interest. In spectral computed tomography, we are interested in the energy channel-wise reconstructions of the object. However, such reconstructions suffer from low signal-to-noise ratio and share the challenges of conventional low-dose computed tomography such as ring artifacts. Ring artifacts arise from errors in the flat-field correction and can significantly degrade the quality of the reconstruction. We propose an extended flat-field model that exploits high correlation in the spectral flat-fields to reduce ring artifacts in the channel-wise reconstructions. The extended model relies on the assumption that the spectral flat-fields can be well-approximated by a low-rank matrix. Our proposed model works directly on the spectral flat-fields and can be combined with any existing reconstruction model, e.g., filtered back projection and iterative methods. The proposed model is validated on a neutron data set. The results show that our method successfully diminishes ring artifacts and improves the quality of the reconstructions. Moreover, the results indicate that our method is robust; it only needs a single spectral flat-field image, whereas existing methods need multiple spectral flat-field images to reach a similar level of ring reduction.


Introduction
Computed Tomography (CT) is a non-invasive imaging technique that allows us to obtain structural knowledge about the interior of objects from a set of projection images. Projection images are acquired by illuminating the object from different angles with radiation from a source e.g., an X-ray beam or beam of neutron radiation. The beam is attenuated as it travels through the object and the attenuated beam is measured by a detector placed opposite the source. The attenuation is governed by absorption in X-ray CT and by scattering in the case of neutron CT. In both cases, the attenuation is material-and energy-specific, and if we measure the attenuation for multiple energies, also referred to as spectral CT, we can obtain a material decomposition of the object.
To obtain the material decomposition, we need to compute the energy channel-wise reconstructions. Most reconstruction methods rely on the assumption that the detector response is known. In practice, however, the detector response is subject to various errors and must be estimated from measurements acquired without an object in the scanner, i.e., from flat-fields, also referred to as air scans [1], white fields [2] or open beams. The flat-fields are noisy due to factors such as measurement noise, miscalibration, defective pixel elements with non-linear response, dust on scintillators, and may introduce concentric rings in the reconstruction, also known as ring artifacts [3]. Ring artifacts are a great challenge for experimental CT set-ups with low-dose and/or short exposure time [4] and can significantly degrade the quality of the reconstruction. In spectral CT, we measure spectral flat-fields, i.e., flat-fields for each energy. However, the spectral measurements share the characteristics of low-dose CT since each energy channel has a low signal-to-noise ratio (SNR) and thus ring artifacts present a challenge in spectral CT [5,6].
To illustrate the challenges of spectral CT, let us consider a neutron CT data set [7] which is described in detail in Section 3. Filtered back projection (FBP) reconstructions of the neutron data are shown for two energies in Figure 1 and the reconstructions reveal presence of ring artifacts. 1.1. Existing methods for ring reduction Several reconstruction methods have been proposed to combat ring artifacts as part of the spectral reconstruction step. Wu et al. [8] propose a reconstruction method that exploits the similarity across spectral images by computing a polychromatic reconstruction (average across the spectral dimension) as a reference image combined with total variation (TV). Lv et al. [9] and Fang et al. [10] both propose deep learning approaches to suppress noise and remove ring artifacts for spectral CT. However, all methods rely on computationally expensive algorithms where the ring reduction is part of the reconstruction estimation process.
Conventional preprocessing methods for ring reduction in monochromatic CT, i.e., single energy CT, can also be applied to the spectral CT data [11,12]. The main drawback here is that the ring reduction techniques must be applied channel-wise to the measured data. Figure 2 shows the eight measured spectral flat-fields stacked vertically and the corresponding singular values. A visual inspection suggests that the spectral flat-fields carry significant redundant information and that we can improve the SNR level in the spectral flat-fields by approximating the spectral flat-fields with a low-rank matrix. In particular, the singular values indicate that the spectral flat-fields can be well-approximated by a rank-one matrix due to the large jump in magnitude between the first and second singular values. A similar idea where principal component analysis (PCA) is used to reduce ring artifacts in case of beam instability has been proposed by Hagemann et al. [13] and Nieuwenhove et al. [14]. The underlying assumptions in these studies are related, but the nature of the problems solved differs. Inspired by Figure 2, we propose an extended flat-field model that exploits high correlation across channels in the spectral flat-fields to reduce ring artifacts in the reconstructions. The extended model relies on the assumption that the spectral flat-fields can be well-approximated by a low-rank matrix. Our method does not depend on a specific reconstruction method since it works directly on the spectral flat-fields. Hence, our method can be combined with all existing reconstruction models such as the conventional FBP or more advanced spectral methods e.g., [8,15,16,17]. Moreover, our method does not need to be applied channel-wise in the sense that the low-rank spectral flat-field simply replaces the measured spectral flat-field in the reconstruction step.

Outline
Section 2 introduces the spectral CT model, existing methods for ring reduction and the proposed methodology. In Section 3, the experimental set-up for the neutron data set is described and numerical experiments are conducted in Section 4. Section 5 discusses the results and Section 6 concludes the paper.

Notation
The set R n is the n-dimensional Euclidean space, R n + is the non-negative orthant and R m×n denotes the set of real-valued m×n matrices. The vector 1 n ∈ R n is a vector of ones, I n×n ∈ R n×n denotes the identity matrix and the transpose of A is denoted A T . The function exp with a vector or matrix as argument are to be interpreted element-wise. If x is a vector, then diag(x) is the diagonal matrix with the elements of x on the diagonal. The 2-norm of a vector is denoted · 2 and the Frobenius norm is denoted · F .

Methods
Consider a spectral data set with m energy channels and let E k denote the kth energy channel. The incident intensity of a beam with energy E k on a detector elemenet is prescribed by the Beer-Lambert law [18], where I(E k ) and I 0 (E k ) are the energy-dependent intensity incident on the detector element and in the object, respectively. Further, is the line segment between the source and the detector, and µ : R d × R → R + is the energy-dependent spatial attenuation function. Let Y k ∈ R rp denote the measurements for the k'th energy with r detector elements and p projection images and discretize the domain into n pixels, then by appropriate discretization (e.g., see [19]), we can describe the measurements by, where ⊗ is the Kronecker product, Z k ∈ R r is energy-dependent intensity incident on the detector, X k ∈ R n are the unknown attenuation coefficients and A ∈ R rp×n is the system matrix describing the traveled distance of the beams. In practice, the energy-dependent intensity incident on the detector is estimated by measuring the detector response. We assume that s spectral flat-fields are measured for each energy. The estimate for Z ∈ R r×m is then conventionally computed by the mean of the s measured spectral flat-fields, i.e., where F = [F T 1 , F T 2 , . . . , F T s ] T and F j ∈ R r×m is the j'th spectral flat-field. However, if F is noisy, then the estimateẐ k can give rise to ring artifacts in the reconstruction.

Low-rank approximation
Essentially, each of the spectral flat-fields carries information about the detector response for all energies and the aim is to exploit the high correlation in the spectral dimension motivated by the observations in Figure 2.
The nearest low-rank matrix of F in the spectral norm can be computed by means of a singular value decomposition (SVD) [20]. An SVD of F is a decomposition where the columns of U ∈ R rs×rs and V ∈ R m×m are orthogonal matrices, and Σ ∈ R rs×m is a matrix with the singular values σ 1 ≥ σ 2 ≥ · · · ≥ 0 on the main diagonal and zeroes elsewhere. The best rank-l approximation of F in the spectral norm is then given by where U i and V i are the ith columns of U and V , respectively. The relative approximation error is given by, which follows from the Eckart-Young-Mirsky theorem.
The rank-one and rank-five approximations of the spectral flat-fields for the neutron data are shown in Figure 3. The approximation error for the rank-one and rank-five matrices are 0.030 and 0.028, respectively. Hence, we only reduce the approximation error with 0.002 by including four extra singular vectors. Considering Figure 3, we see that the rank-one approximation has mitigated a substantial amount of noise whereas in the rank-five approximation, some of the noisy tendencies in the spectral flat-fields start to reappear. This observation is confirmed by considering the difference images between the spectral flat-fields and the low-rank matrices in Figure 3e and 3f. Hence, for the numerical experiments, we will confine ourselves to considering the rank-one approximation of the spectral flat-fields, i.e., we compute the estimate forẐ by replacing F by At this point, we emphasize that our method is computationally cheap as it scales with O(min(rs, m) 2 max(rs, m)). Moreover, it can be combined with any reconstruction method, e.g., FBP, iterative methods, statistical models, etc. since it is only applied to the spectral flat-fields. The corresponding difference images between the spectral flat-fields and low-rank matrices are shown in (e) and (f) for rank-one and rank-five, respectively.

Existing ring reduction method
Our method works solely on the spectral flat-fields and thus separates the ring reduction from the reconstruction step. Hence, we compare our method to two existing ring reduction techniques for monochromatic CT. We compare our methodology with the preprocessing method proposed by Münch et al. [11] which combines wavelet and Fourier filtering to mitigate ring artifacts in the reconstruction. The method is computationally inexpensive and does not increase the overall computational cost significantly. The method depends on three parameters; we use a damping factor of 0.9 and the Daubechies 5 wavelet with a three-level decomposition for all numerical experiments, see [11] for details. The second preprocessing method is proposed by Vo et al. [12] and uses a combination of sorting and smoothing (non-local means) in attempt to smooth the data and thereby reduce ring artifacts. The method depends on a parameter related to the smoothing filter, we choose a parameter value of 31, see [12] for further details. We denote the methods from Münch et al. [11] and Vo et al. [12] by FW (Wavelet Fourier) and NLM (non-local means), respectively.

Neutron data
We validate the proposed methodology on a neutron CT data set [7]. The neutron data were acquired at the imaging and materials science and engineering (IMAT) beamline operating at the ISIS spallation neutron source (Rutherford Appleton Laboratory, UK). Figure 4 shows a sketch of the object of interest. The object consists of six cylinders whereof five are filled with metal powders, i.e., aluminum (Al), iron (Fe), copper (Cu), nickel (Ni) and zinc (Zn) powders.

Zn
Ni Fe Al Cu Figure 4. Sketch of the object for the experimental neutron data [7]. The experimental set up consists of six cylinders where five of the cylinders have been filled with aluminum (Al), iron (Fe), copper (Cu), nickel (Ni) and zinc (Zn) powders. The sixth cylinder is empty.

Data acquisition
Ametova et al. [17] describe the data acquisition and pre-treatment of data in detail, i.e., beam instabilities, overlap correction and spectral averaging. We consider the reconstruction-ready measurements and confine ourselves to reporting only the essential details of the experimental set-up. The data set contains m = 339 spectral projections acquired at p = 120 equidistant angles distributed from 0 • to 180 • with 1.5 • angular increments. Eight spectral flat-fields were acquired, four prior to the scan and four after, i.e., s = 8. Each projection consists of 460 × 460 pixels with a pixel size of 0.055 mm resulting in a view of approximately 25 × 25 mm 2 and r = 460 detector elements and n = 460 2 pixels for the experiments. We choose the 127th vertical detector row for experiments, i.e., we consider a two-dimensional set-up. The experimental setup is summarized in Table 1. The attenuation coefficient for the neutron experiment is denoted Σ tot (λ) and has unit cm −1 , see [17] for further details.

Numerical experiments
We consider two reconstruction models: FBP and a weighted least squares (WLS) reconstruction model combined with TV regularization. We compare our method to the conventional flat-field correction and the existing ring reduction techniques described in Subsection 2.2. Table 2 provides an overview of the reconstruction models and ring reduction techniques. WF-FBP WF-TV Preprocessing Vo et al. [12] NLM-FBP NLM-TV Low-rank spectral flat-fields LR-FBP LR-TV We use AIR TOOLS II [21] to generate the parallel-beam geometry of the experimental setup. For the FBP reconstructions, we use the fbp function from AIR TOOLS II with the Hann filter to reduce the noise in the computed FBP reconstructions. For the TV reconstructions, we used the implementation of WLS with TV from [4] in MATLAB. All TV reconstructions have a maximum number of iterations of 1000 and a regularization parameter of 0.005. The regularization parameter was found by visual inspection of reconstructions for varying values of the regularization parameter.
Note that we pick the same regularization parameter for all energies which will result in some reconstructions being a bit over-regularized whereas other reconstructions might be slightly under-regularized. The reason is that the SNR changes significantly as a function of energy and thus there is not a single regularization parameter that fits all energy channels. However, the purpose of the methods is ring reduction and thus we limit the experiments to considering the same regularization parameter for all reconstructions.

Error measures
The quality of the computed reconstructions is assessed by visual inspection combined with contrast-to-noise ratio (CNR). The CNR metric is used for evaluating the image contrast and noise properties for a selected region of interest (ROI) [22]. We use the method proposed by Bian et al. [23] where a ROI with low-contrast structure is compared to a background ROI while taking the standard deviations of both the signal and background ROIs into account. The ROIs are clearly marked on the figures when applicable.

Experiment: Different energies
There are too many energies to visualize all reconstructions, and the SNR varies significantly between the energies [17]. Thus, we select three energies based on the relative difference (RD) between the computed FBP and LR-FBP solution to ensure a representative visualization of energy channels. We define RD by the measure, where X FBP k and X LR-FBP k denote the FBP and LR-FBP reconstructions for energy k, respectively. We select three energies corresponding to the minimum, median and maximum RD, corresponding to energies 3.1Å, 4.2Å and 2.9Å, respectively. The FBP reconstructions in Figure 5 (first column) are all affected by ring artifacts. However, the severity of the ring artifacts increases from top to bottom, i.e., when the RD increases. The WF-FBP reconstructions (second column) seem to mitigate the ring artifacts at first glance, however, when carefully comparing the WF-FBP with NLM-FBP and LR-FBP, wavelike artifacts can be seen in the WF-FBP reconstructions. Hence, the ring artifacts have been reduced but not eliminated. For the NLM-FBP and LR-FBP reconstructions (third and fourth column), there are no visible artifacts remaining.  Considering the CNR listed in Table 3, we see that the CNR decreases from top to bottom for all four methods. Thus, when the SNR is low, our method differs most from the FBP, which can be explained by the fact that the ring artifacts are more dominating when the SNR is low. The CNR for all four methods is quite close in Figure 5. This might be explained by the fact that the dominating noise contribution comes from the measurement noise and not the ring artifacts, since all FBP reconstructions suffer from a very low SNR. Note that the high noise level might conceal remaining ring artifacts.
Inspecting the TV reconstructions in Figure 6, we see that TV regularization reduces the noise level significantly, and consequently, the ring artifacts appear more severe for the TV reconstructions. In addition, we also note that the vague wave structures in the WF-FBP reconstructions are even clearer in the WF-TV reconstructions compared to the WF-FBP reconstructions. By carefully inspecting the NLM-TV reconstruction for 2.9Å, one can see that the preprocessing has introduced a dark spot with negative values in the center. The other NLM-TV reconstructions show no introduced artifacts and closely resemble the LR-TV reconstructions (fourth column). The LR-TV reconstructions reveal no ring structure even though the noise level is reduced.
The TV reconstructions have significantly improved the overall CNR compared to the FBP reconstructions, which can be seen in Table 3. For the TV experiment in Figure 6, we also see an increase in CNR when applying a ring reduction method, especially for NLM-TV and LR-TV. The LR-TV reconstruction achieves the highest CNR for two out of the three energies depicted in Figure 6.

Experiment: Effect of number of spectral flat-fields
We now perform an experiment with varying number of spectral flat-fields to validate the robustness of our proposed method. We consider energy 2.1Å for the experiment and use eight, four, two and one flat-field for the flat-field correction, respectively. The FBP, WF-FBP, NLM-FBP and LR-FBP reconstructions are depicted in Figure 7 and the TV reconstructions for the same experiment are shown in Figure 8. The CNR measures for the reconstructions are listed in Table 4.
If we start by considering the first column of Figure 7, we see that there are visible ring artifacts in all FBP reconstructions, but the severity of the ring artifacts increases as the number of flat-fields used for flat-field correction decreases (i.e., from top to bottom). The WF-FBP Figure 8. TV, WF-TV, NLM-TV and LR-TV reconstructions of energy 2.1Å with eight, four, two and one flat-field, respectively. The white squares with full and dashed lines mark the structure and background ROIs, respectively. Table 4. CNR for experiment depicted in Figure 7 and Figure 8  reconstructions are shown in the second column of the figure. Vague ring artifacts can be seen in the WF-FBP reconstruction using all eight flat-fields. However, the severity of these ring artifacts increases when the number of flat-fields decreases, just as for the FBP reconstructions.
The NLM-FBP reconstructions in third column show less wave-like artifacts than NLM-FBP, however, then using four or fewer flat-fields, the method starts to struggle, and artifacts start to arise in the reconstructions. The LR reconstructions in the fourth column show no sign of artifacts, not even in the case where a single flat-field is used and thus our method seems to be robust. The robustness of LR-FBP can also be seen in the CNR reported in Table 4. Here the CNR is almost constant for the four reconstructions with LR, whereas the CNR increases with the number of flat-fields for both FBP, WF-FBP and NLM-FBP. The experiment was repeated using TV for the reconstruction model and the findings support those found using FBP. The SNR is significantly increased by using TV regularization and the remaining artifacts for the preprocessing methods are even more pronounced in Figure 8, which is also supported by the CNR measure listed in Table 4.

Discussion
An explanation for the robustness of our method lies with the fact that the spectral dimension carries lots of redundant information. For the experiment, we have eight flat-fields but 339 energy channels. Thus, using a single spectral flat-field corresponds to 339 measurements for our method, whereas eight flat-fields with the conventional method only corresponds to eight flat-fields. Thus, if the assumption of low rank holds, then more energy channels will give more redundancy in the spectral flat-fields and thereby our proposed method will gain further advantage over the existing methods. We investigated spectral plots for material decomposition, i.e., plot of a single pixel across the spectral dimension for the distinct materials. However, preliminary results showed no visual effect of applying the preprocessing methods.
We could consider the spectral flat-fields as a three-dimensional tensor with dimensions s × r × m. While tensor decomposition methods exist, such as the Tucker decomposition and parallel factors decomposition (PARAFAC) [24,25,26], preliminary experiments showed no significant difference between the low-rank tensor, and the low-rank matrix obtained by treating the spectral flat-fields as a matrix as in (3).
A possible improvement of the model could be the inclusion of the low-rank matrix as parameters of the reconstruction model and jointly estimating the reconstruction and the spectral detector response, e.g., see [19]. A natural initial guess for the low-rank matrix would be the estimate obtained by the LR method.

Conclusion
We have proposed an extended flat-field model for spectral CT that exploits high correlation in the spectral flat-fields by replacing the spectral flat-fields with a low-rank approximation to mitigate ring artifacts. The proposed methodology can be combined with any existing reconstruction method and only depends on a single parameter, which is easy to choose by inspection of the singular values of the spectral flat-fields. The method was validated on neutron CT data and compared to the conventional flat-field correction and two preprocessing methods for ring reduction. Our method successfully mitigated ring artifacts in all experiments whereas the other methods struggled to suppress ring artifacts, especially in more challenging cases with severe ring artifacts.