Analytical expressions for the reconstructed image of a homogeneous cylindrical sample exhibiting a beam hardening artifact in X-ray computed tomography

Abstract

BACKGROUND:

Cylindrical phantoms are often imaged by X-ray computed tomography (CT) to evaluate the extent of beam hardening (or cupping artifact) resulting from a polychromatic X-ray source.

OBJECTIVE:

Our goal was to derive analytical expressions for the reconstructed image of a homogeneous cylindrical phantom exhibiting a cupping artifact, to permit a quantitative comparison with experimental cupping data.

METHODS:

A filtered backprojection method was employed to obtain the analytical cupping profile for the phantom, assuming that the projection data could be approximated as a power series with respect to the sample penetration thickness.

RESULTS:

The cupping profile was obtained analytically as a series of functions by employing Ramachandran filtering with an infinite Nyquist wavenumber. The quantitative relationship between the power series of the projection and the nth moment of the linear attenuation coefficient spectrum of the phantom was also determined. Application of the obtained cupping profile to the evaluation of the practical reconstruction filters with a finite Nyquist wavenumber and to the best choice of the contrast agent was demonstrated.

CONCLUSIONS:

The set of exact solutions derived in this work should be applicable to the analysis of cylindrical phantom experiments intended to evaluate CT systems.

Keywords

Beam hardening cupping artifact exact solution phantom polychromatic X-ray spectrum X-ray CT

1 Introduction

Homogeneous cylindrical samples are often used as phantoms in X-ray computed tomography (CT) experiments with a polychromatic X-ray source [1–13 , 35]. The resulting CT images containing undesirable cupping artifacts [1–3 , 14] can be used to assess the extent of beam hardening by the CT system. To do so quantitatively, the theoretical line profile of the CT value along the radial direction of the reconstructed two-dimensional (2-D) cylindrical sample image [1 , 15] is required to compare with the experimentally obtained line profile. The theoretical line profile contains undesirable errors due to, for example, the imperfect filtering and discretization [4 , 37] if the profile is obtained by the numerical image-reconstruction simulation. The ideal beam-hardening assessment without such undesirable errors is possible if the analytical (not numerical) line profile is available. However, to the best of our knowledge, there is no paper on the general solution for the analytical line profile.

In the present theoretical study, we therefore derived analytical expressions for the reconstructed image of a homogeneous cylindrical phantom exhibiting a cupping artifact by assuming the following; (i) that non-linear projection data can be considered as a power series with respect to the sample penetration thickness and (ii) that Ramachandran filtering with an infinite Nyquist wavenumber can be employed. We also examined the relationship between the power series for the projection data and the nth moment of the linear attenuation coefficient (LAC) spectrum.

The obtained analytical solutions were applied to the evaluation of the practical reconstruction filters. Numerical simulations reproducing cupping artifacts for a finite Nyquist wavenumber in association with Ramachandran, Shepp, and Chesler filters were performed for a cylindrical specimen containing an aqueous KI solution to confirm the applicability of the analytically-obtained solutions. The analytical solutions were also applied to the choice of the best contrast agent for a cylindrical aqueous phantom. We searched for the atomic number of the heavy element in the aqueous solution/suspension phantom that induced the least cupping artifact.

2 Analytical expressions for a cupping artifact

2.1 From the projection to the cupping profile

Herein, we consider the analytical 2-D image reconstruction of a homogeneous cylindrical sample with radius, R, using a polychromatic X-ray source and an array of linear detectors (Fig. 1). The bowtie filtration [35] which modulates the original X-ray source spectrum is not considered. The incident parallel X-rays have an initial intensity, I(0), and penetrate the sample to a depth of s. The attenuated X-ray intensity, I(s), is measured physically [36] by detectors placed at a constant distance of δ. The projection, p(s), obtained by the detectors is defined as:

Fig. 1

The configuration of a cylindrical sample (in grey) with radius R in an imaging system with a detector gap δ. The polar coordinates of an arbitrary point (the yellow dot) within the sample are (r, θ). Another coordinate system (an x - y system) is rotated by 180 degree during the X-ray emissions. The sample center coincides with the center of rotation.

p (s) = ln [\frac{I (0)}{I (s)}] for 0 \leq s \leq 2 R .

(1)

The 2-D polar and Cartesian coordinate systems, (r, θ) and (x, y), respectively, are indicated in Fig. 1. The system is axisymmetric because the sample is homogeneous and cylindrical. Thus, 2-D image reconstruction by CT imagery is essentially identical to deriving the 1-D cupping profile, f (r), rather than f (x, y) where f is the reconstructed value of the LAC obtained by the CT system.

Typical sample LAC and X-ray source spectra are provided in Fig. 2a. The energy-dependent nature of the spectra in Fig. 2a results in a non-linear projection, p(s), and non-uniform reconstructed value, f(r), within the phantom. With regard to the non-linearity of the projection, we assume that p(s) can be suitably approximated by h(s), which is a power series with respect to the penetration thickness, s, written as [1]:

Fig. 2

A comparison between the reconstructed voxel values obtained from the analytical expression and a simulated CT image with respect to a homogeneous cylindrical phantom (a 440 mM aqueous KI solution). (a) The normalized photon-energy spectrum of the X-ray source (at an acceleration voltage of 100 kV). Two LAC spectra of the KI solution and of a 160 mM samarium aqueous solution/suspension are superimposed. The K absorption edges of iodine and samarium are indicated by arrows (33 and 47 keV, respectively). (b) The reconstructed CT image of the cylindrical KI solution phantom obtained using the Ramachandran filter with a finite Nyquist wavenumber. The image dimensions were trimmed to 300² voxels (= 3² cm²), and the radius, R, of the phantom is 90 voxels (= 0.9 cm). The yellow line (the r-axis) expanding from the sample center to the image edge is the baseline for the line profile. (c) Line profiles of the reconstructed voxel values along the yellow line in (b) for three reconstruction filters (Ramachandran, Shepp, and Chesler). The true line profile calculated using the exact solution based on Equations (3) and (4) is also plotted for N = 100. (d) An expanded view of (c) around the sample edge at r = 0.9 cm. (e) The projection, p, as a function of s as calculated using Equations (1) and (7). The projection, h, approximated as a power series of s as calculated using Equation (2) via μ _n (0) from Equation (A9) is also shown. The linear component of Equation (2) (i.e., C₁ s) is superimposed by a dotted line to depict the non-linearity of p and h.

h (s) = \sum_{n = 1}^{N} C_{n} s^{n} for 0 \leq s \leq 2 R,

(2) where C_n is a coefficient and N is an integer representing the upper limit of the summation. The energy-dependent LAC spectrum and polychromatic X-ray source spectrum yield the cupping artifact. The physical information of the two spectra is combined and embedded in Equation (2). Thus, it is possible to derive the analytical expressions of the cupping profile by the conventional image-reconstruction algorithm using the non-linear projection data of Equation (2).

Employing the conventional filtered backprojection method [4 , 17], we derived an analytical expression for the cupping profile, f(r), as a series of functions, written as: $f (r) = \sum_{n = 1}^{N} {F_{n} (R^{2} - r^{2})}^{\frac{n - 1}{2}} for 0 \leq r < R$ (3)

and $f (R) = \frac{C_{1}}{2} .$ (4)

Here, the coefficient F_n is given by: $F_{n} = \frac{2^{n}}{\sqrt{π}} \frac{Γ (\frac{n}{2} + 1)}{Γ (\frac{n + 1}{2})} C_{n},$ (5)

where Γ is the gamma function. The detailed derivation of Equations (3) and (4) from Equation (2) is described in Appendix 1. It should be noted that Ramachandran filtering with an infinite Nyquist wavenumber (i.e., infinitesimal δ) was employed in this derivation.

Equation (3) shows that the linearity (or non-linearity) of the projection is transferred to the constant (or r-dependent) distribution of f(r) as follows. (i) The coefficient for the linear component in Equation (2), C₁, is embedded as F₁ in the right-hand side of Equation (3), which contributes to the constant component of f(r) without cupping artifact. (ii) Each coefficient C_n for n≥2 corresponding to the non-linear term in Equation (2) occurs in F_n on the right-hand side of Equation (3) as the r-dependent component of f(r).

2.2 From the spectral moment to the projection

The non-linearity of the projection in Equation (2) is a consequence of the energy-dependence of the LAC spectrum and the broadness of the X-ray source spectrum. We therefore developed a quantitative relationship between the projection and spectra by introducing a spectral moment. The incident X-ray intensity, I(0), in Equation (1) can be considered as the sum of all components having various photon energy values, E, and thus can be summarized as: $I (0) = \int_{0}^{\infty} i_{0} (E) dE,$ (6) where i₀(E) is the energy-dependent spectrum of the X-ray source. Hereafter, we assume that I(0) is normalized to a value of unity. Following the penetration of the cylindrical sample in Fig. 1 by the X-ray beam, the X-ray intensity is attenuated by Beer’s law to [18]: $I (s) = \int_{0}^{\infty} i_{0} (E) exp (- μ (E) s) dE,$ (7) where μ(E) is the energy-dependent LAC spectrum of the cylindrical sample.

In order to relate the power-low non-linearity in Equation (2) to the exponential non-linearity in Equation (7), we introduce the nth moment of the LAC spectrum, weighted by the attenuation due to the sample penetration, defined as: $M_{n} (s) = \int_{0}^{\infty} i_{0} (E) μ (E)^{n} exp (- μ (E) s) dE for n = 0, 1, 2, \dots, N,$ (8) where M_n (s) is the nth moment. M_n (s) can be normalized by I(s) as: $μ_{n} (s) = \frac{M_{n} (s)}{I (s)} for n = 0, 1, 2, \dots, N .$ (9)

Through the process described in Appendix 2, μ_n (0) can be related to C_n as the following recurrence relation: $C_{n + 1} = - ν_{n + 1} - \sum_{m = 1}^{n} ν_{n - m + 1} \frac{m}{n + 1} C_{m} for n = 0, 1, 2, \dots, N - 1,$ (10)

where the coefficient ν _n is defined by: $ν_{n} = \frac{{(- 1)}^{n} μ_{n} (0)}{n!} = \frac{{(- 1)}^{n}}{n!} \int_{0}^{\infty} i_{0} (E) μ (E)^{n} dE for n = 0, 1, 2, \dots N .$ (11)

By the above process, we obtained the analytical relationships among the normalized spectral moment (μ _n (0) in Equation (A9)), the coefficient for the non-linear projection (C_n in Equation (2)), and the coefficient for the cupping profile (F_n in Equation (3)). In the following sections, these exact solutions were applied to the comparison with the practical reconstruction filters and to the search for the best contrast agent in terms of the reduction of the cupping artifact.

3 Application of the exact solutions

3.1 Effects of reconstruction filters

The analytical expressions derived in Sec. 2 assume Ramachandran filtering with an infinite Nyquist wavenumber. This implies that the detector interval, δ, is zero. In practical CT imagery, however, the quantity δ is not zero, and Shepp (i.e., sinc) and Chesler (i.e., hanning) reconstruction filters with finite Nyquist wavenumbers are often used in addition to the Ramachandran filter [4, 16]. Thus, conventional numerical simulations to reconstruct a synthetic CT image exhibiting a cupping artifact were performed using such practical filters to demonstrate the applicability of these analytical expressions to practical CT imagery.

The cylindrical sample selected for the reconstruction simulations was composed of a homogeneous aqueous 440 mM KI solution. The radius of the cylinder, R, was 0.9 cm and the detector separation distance, δ, was 0.01 cm. Thus, the radius of the reconstructed cylindrical sample was (0.9 cm/0.01 cm) = 90 voxels. The LAC spectrum of the sample as calculated using the NIST database [19] is shown in Fig. 2a, in which the X-ray source spectrum at an acceleration voltage of 100 kV [14] is superimposed. The target of the X-ray tube is made of Mo-W alloy. Metal films (Al film of approximately 1 mm and Cu film of 0.1 mm in thickness) are preinstalled on the CTsystem [14].

The 2-D image reconstruction of the cylindrical sample in Fig. 1 was performed with our original computer program. The details of the computer program have been described elsewhere [4 , 20]. Briefly, the parallel X-ray beams with the polychromatic X-ray source spectrum are assumed, and the attenuated X-ray intensity is measured on the 1-D array of detectors. The non-linear projection, p(s), is calculated on the discrete points (separation distance, δ) using Equations (1) and (7) for 0≤s≤2 R. Totally 805 projection data were acquired on the 1-D detector array during the sample rotation of 180°. The original dimensions of the 2-D reconstructed image were 512² voxels, and the obtained sinogram consisted of 512×805 voxels. The convolution backprojection method was applied to the sinogram to reconstruct a 512²-voxel image with voxel dimensions of 0.01² cm². Three typical reconstruction filters, the Ramachandran, Shepp, and Chesler filters, with a finite Nyquist wavenumber of 1/(2δ) = 50 1/cm, were employed in the computer program. The grey scale of the reconstructed TIFF image was 16 bits. We applied a line-profile analysis of the voxel values to the 16-bit image exhibiting beam hardening to evaluate the effects of the reconstruction filters on the degree of cupping artifact. The ImageJ program developed at the National Institutes of Health was used for the line profile analysis of the 16-bit TIFF image.

The exact solutions corresponding to the KI phantom of R = 0.9 cm were calculated as follows. The upper limit of n (i.e., N), was taken to be 100. Using the X-ray source and LAC spectra in Fig. 2a, values for μ_n(0) were calculated via Equation (A9) for n = 1 to 100. These values were subsequently converted to ν_n using Equation (11), and C_n was calculated by Equation (10). Finally, F_n was determined by Equation (5) to obtain the cupping profile, f(r).

3.2 Search for the best contrast agent

Substance with a high atomic number (e.g., iodine and lanthanoids) is used as a contrast agent in X-ray CT [6]. The use of such heavy elements yields the undesirable cupping artifact [9], which makes the quantitative CT image analysis difficult. The degree of the cupping depends on the atomic number of the heavy element [10]. In this section, we applied the exact solutions obtained in Sec. 2 to the search for the atomic number of the best contrast agent that induces the least cupping artifact.

The degree of the cupping artifact of a cylindrical phantom can be expressed as the difference in the f (r) values between the center and the rim of the reconstructed phantom image. As for the f (r) value at the rim, the limiting behavior of Equation (3) described in Appendix 1 ensures the approximation that f (r) ≈C₁ in the immediate vicinity of the sample rim. Thus, we take that the best contrast agent yields a minimum with respect to the C₁ value for a fixed f (0) value.

We consider the following simple case. A mixture of simple substance of a heavy element and pure water is supposed as a cylindrical homogeneous phantom. The mixture is an aqueous solution or a fine-grained suspension, so that the grain size of the simple substance is very small compared with the voxel dimension, and the mixture can be well approximated as homogeneous.

Heavy elements ranging from ₄₉In to ₇₉Au were examined. The f(0) value was tentatively fixed to be 0.48 1/cm (≈ 1000 Hounsfield unit (HU)) in the present study. The detailed input data for the exact solutions were as follows. The radius of the phantom, R, was 0.9 cm, N= 100, and the X-ray source spectrum was for the 100 kV acceleration, which was the same as the simulation in Sec. 3.1. We obeyed Nakashima and Nakano [10, 20] for the estimation of the bulk density of the mixture, which was needed for the LAC spectrum calculation. An example for the calculated LAC spectrum is shown in Fig. 2a for the samarium-water mixture system.

4 Results

4.1 Effects of reconstruction filters

The synthetic 2-D CT image of the cylindrical phantom of the 440 mM KI solution is shown in Fig. 2b. The original 512² voxels image was trimmed to 300² voxels to provide an expanded view of the cylindrical sample. This figure demonstrates a characteristic beam hardening artifact (cupping phenomenon) along with a non-uniform distribution of the reconstructed LAC values within the sample (see the bright voxels near the sample rim). The line profile along the radial direction (the yellow baseline in Fig. 2b) is provided in Fig. 2c, and also shows the cupping phenomenon. Figure 2c is enlarged in Fig. 2d to enhance the differences between the three reconstruction filters.

With regard to the analytical solutions, the theoretical cupping profile, f(r), calculated using Equations (3) and (4) is presented in Fig. 2 c and d. A data point spacing of 0.001 cm was applied in the case of this theoretical plot, which is an order of magnitude smaller than the detector spacing (i.e., & δ = 0.01 cm) used for the simulation. This was done because the theoretical profile is for an infinite Nyquist wavenumber (an infinitesimal & δ). The μ_n(0), ν_n, C_n, and F_n values calculated for Fig. 2 are summarized in Table 1 for the n values up to ten. The two non-linear projections (p(s), directly calculated by Equations (1) and (7), and h(s), indirectly calculated by Equation (2) via μ_n(0)) are plotted in Fig. 2e, and demonstrate reasonable agreement.

Table 1
List of the coefficients, μ _n (0), ν _n , C_n, and F_n in Fig. 2 for the n values up to ten. The units of these four coefficients are common, but depend on n (i.e., 1/cm ⁿ ). The two values for n = 0 (i.e., μ₀(0) = 1 and ν₀ = 1) are dimensionless constants that are independent of the size and chemistry of the sample

n μ _n (0) ν _n C_n F_n

0 1 1 N/A N/A

1 0.96208 −0.96208 0.96208 0.96208

2 1.14125 0.57062 −0.10783 −0.27458

3 1.60713 −0.26786 0.01570 0.09421

4 2.56714 0.10696 −0.00045 −0.00605

5 4.47574 −0.03730 −0.00056 −0.01666

6 8.28798 0.01151 0.00014 0.00920

7 16.01007 −0.00318 2.4 e-7 3.4 e-5

8 31.88811 0.00079 −8.3 e-6 −0.00247

9 64.98430 −0.00018 1.8 e-6 0.00113

10 134.79017 3.7 e-5 1.2 e-7 0.00016

n	μ _n (0)	ν _n	C_n	F_n
0	1	1	N/A	N/A
1	0.96208	−0.96208	0.96208	0.96208
2	1.14125	0.57062	−0.10783	−0.27458
3	1.60713	−0.26786	0.01570	0.09421
4	2.56714	0.10696	−0.00045	−0.00605
5	4.47574	−0.03730	−0.00056	−0.01666
6	8.28798	0.01151	0.00014	0.00920
7	16.01007	−0.00318	2.4 e-7	3.4 e-5
8	31.88811	0.00079	−8.3 e-6	−0.00247
9	64.98430	−0.00018	1.8 e-6	0.00113
10	134.79017	3.7 e-5	1.2 e-7	0.00016

The effects of the N value in Equation (3) for the case of Fig. 2 are shown in Fig. 3. The f(r) value is constant (i.e., f(r) = C₁ = μ₁(0)) and no beam hardening occurs if N = 1, which is a consequence of assuming that the projection, h(s), in Equation (2) is linear with respect to s. From this figure it is evident that Equation (3) converges rapidly on a fixed profile in conjunction with a small degree of oscillation as N increases.

Fig. 3

The dependence of the exact solution, f(r), on N. Here, f(r) for N = 100 is the same as that plotted in Fig. 2c, and it is difficult to distinguish the f(r) curves for N = 6, 10, and 100 due to the good agreement.

4.2 Search for the best contrast agent

Using the X-ray source spectrum in Fig. 2a and a LAC spectrum of a mixture of a specific heavy element and water, values for μ_n(0) were calculated via Equation (A9) for n = 1 to 100. These values were subsequently converted to ν_n using Equation (11), and C_n was calculated by Equation (10). Finally, F_n was determined by Equation (5) to obtain the cupping profile, f (r). By trial and error, we searched for the molar concentration of each heavy element in the aqueous solution/suspension which yields f (0)= 0.48 1/cm (≈ 1000 HU), and plotted the C₁ value against the atomic number of the heavy element as well as the corresponding molar concentration. The results are summarized in Fig. 4, depicting that the C₁ value is sensitive to the atomic number of the heavy element in the mixture. For example, C₁= 0.522 (≈ 1180 HU), 0.500 (≈ 1080 HU), and 0.545 (≈ 1270 HU) for ₅₃I, ₆₂Sm, and ₇₉Au, respectively. This is a consequence of that the degree of the cupping artifact is sensitive to the relative position between the K absorption edge of the heavy element and the spectrum peak of the X-ray source [9 , 20].

Fig. 4

Dependence of C₁ value on the atomic number of the heavy element suspended in the water sample for a fixed f (0) value of 0.48 1/cm. The corresponding molar concentration of the heavy element is also plotted. Note that samarium (atomic number, 62) yields the minimum cupping artifact of C₁ – f (0) ≈ 0.02 1/cm, for which the corresponding LAC spectrum is shown in Fig. 2a.

5 Discussion

From Fig. 2c, it can be seen that the line profile of each simulated image closely coincides with the exact solution. This result confirms that the power-series approximation for the non-linear projection (Equation 2) is reasonable. This is supported by Fig. 2e, depicting good agreement between h(s) and p(s). As for the choice of the upper limit of the n value (i.e., N), a larger N value theoretically ensures better accuracy as far as the numerical calculation error is not serious. However, the three profiles for N = 6, 10, and 100 yield almost the same curves in Fig. 3, implying that an N value as small as approximately 6 would be practically sufficient for the KI-solution sample of Fig. 2. This is also confirmed by the rapid decrease of the F_n values with increasing n in Table 1.

Figure 2d shows a discrepancy in the line profiles of the simulated images and the exact solution near the sample rim (i.e., at r = 0.9 cm). This is attributed to the use of a finite Nyquist wavenumber when generating the simulated image [4, 16]. The discrepancy relative to the exact solution is greatest for the Chesler filter and smallest for the Ramachandran filter. The Chesler and Ramachandran filters are the strongest and weakest low-pass filters among the three options, respectively, and this is reflected in the degree of discrepancy in Fig. 2d. However, the three simulated image line profiles are in reasonable agreement with the exact solution, with the exception of a few voxels at the sample rim (Fig. 2c and d). Thus, we conclude that the exact solution is applicable to the analysis of real-world CT images using either Chesler, Shepp, or Ramachandran filters with finite Nyquist wavenumbers, except in the vicinity of the sample rim.

Figure 4 depicts that the best contrast agent is not conventionally used iodine [9] but samarium [6] because the difference of the f (r) values between the center and the rim is as small as 0.020 1/cm (≈ 80 HU) for samarium while it is 0.042 1/cm (≈ 180 HU) for iodine. This is a consequence of that the K absorption edge of samarium (47 keV) is located near the spectrum peak of the continuous X-ray component of the source spectrum (see Fig. 2a) [9 , 20]. It should be noted that Fig. 4 was obtained not by the conventional image-reconstruction numerical simulations having undesirable errors [4 , 37] but by the analytical solutions. The analytical solutions derived in the present study would contribute to the reliable assessment of the contrast agent performance without such undesirable errors.

As can be seen from Equations (11), (10), and (5), the nth moment of the LAC (i.e., μ_n(0)) and the non-linearity of the projection (i.e., C_n) are mathematically related to the coefficient, F_n, for the line profile, f(r). This set of analytical expressions can be considered as either a forward or inverse problem. Assuming a forward problem [4 , 20], it is possible to predict the line profile, f(r), for the reconstructed CT image of a sample with a known size and chemistry if the X-ray source spectrum or thickness-dependent projection data are available. This would have applications in assessing the performance of a CT imaging system using a cylindrical phantom as demonstrated in Fig. 4. In the case of an inverse problem [18 , 22], the same set of analytical expressions allows us to estimate the X-ray source spectrum by analyzing the line profile data, f(r), experimentally obtained using cylindrical phantoms of known size and chemistry. In term of the mathematical inverse problem, the use of plural phantoms with different heavy elements is preferable because it would increase the estimation accuracy of the X-ray source spectrum. For example, the C₁ values for ₅₈Ce and ₇₄W are significantly different in Fig. 4, and both can be independent constraint in solving the inverse problem. Thus, the simultaneous use of the ₅₈Ce- [20] and ₇₄W-water [9] phantoms would contribute to the more accurate estimation of the source spectrum.

The exact solutions obtained are applicable not only to X-ray CT [13 , 34], but also to other diverse imaging systems employing straight ray paths. Example includes those employing neutrons, positrons, muons, and neutrinos [27 –30]. The analytical expressions for a homogeneous cylindrical phantom could also be used to assess imaging systems for a wide range of CT applications.

6 Conclusions

Analytical expressions were obtained for the reconstructed image of a homogeneous cylindrical sample exhibiting beam hardening (or cupping artifact) due to a polychromatic X-ray source. The quantitative relationship between the cupping profile and the LAC and X-ray source spectra was derived from non-linear projection data with respect to the penetration thickness by employing Ramachandran filtering with an infinite Nyquist wavenumber. The cupping profile resulting from the analytical expressions was in reasonable agreement with that generated by numerically simulated CT imagery in conjunction with various reconstruction filters (Ramachandran, Shepp, and Chesler filters with a finite Nyquist wavenumber), except for a few voxels near the sample rim. Application of the analytical expressions to the best choice of the contrast agent that induced the least cupping artifact was also demonstrated successfully. The exact solutions as determined in the present study are expected to have applications in the analysis of cylindrical phantom data during the evaluation of real-world CT systems using such practical reconstruction filters.

Conflict of interest

The authors declare that they have no conflict of interests.

Appendix 1: Filtered backprojection for a homogeneous cylindrical sample

The derivation of Equations (3) and (4) from Equation (2) using the conventional filtered backprojection method [4 , 17] is described herein. First the one-dimensional Fourier transform of the projection was calculated, second the Ramachandran filtering followed by the inverse Fourier transform was performed, and finally the backprojection operation was employed for the reconstruction of the one-dimensional image, f(r).

Applying a one-dimensional Fourier transform to the projection h(s), the quantity P(z) is obtained as: $\begin{matrix} P (z) = \int_{- \infty}^{\infty} h (s) exp (- 2 π izt) dt \\ = 2 \int_{0}^{R} h (2 \sqrt{R^{2} - t^{2}}) cos (2 π zt) dt, \end{matrix}$ (A1) where z is the wavenumber and the parameter t (= rsinθ) is as indicated in Fig. 1. The evenness of h with respect to t is also incorporated into the above equation.

Equation (2) was substituted into Equation (A1) and the integral formula (ET I 11(8)) of [31] was used to obtain: $\begin{matrix} P (z) = \sum_{n = 1}^{N} C_{n} 2^{n + 1} \int_{0}^{R} {(R^{2} - t^{2})}^{\frac{n}{2}} cos (2 π zt) dt \\ = \sum_{n = 1}^{N} C_{n} {(\frac{2}{\sqrt{π}})}^{n} {(\frac{R}{z})}^{\frac{n + 1}{2}} Γ (\frac{n}{2} + 1) J_{\frac{n + 1}{2}} (2 π Rz), \end{matrix}$ (A2)

where J_(n+1)/2 (2πRz) is the Bessel function of the first kind.

Ramachandran filtering (i.e., |z|) with an infinite Nyquist wavenumber was applied to P(z) to obtain the quantity q(t), written as: $\begin{matrix} q (t) = \int_{- \infty}^{\infty} P (z) | z | exp (2 π itz) dz \\ = 2 \int_{0}^{\infty} P (z) z cos (2 π tz) dz . \end{matrix}$ (A3)

The evenness of the integrand was incorporated into the above equation.

As a final step, the backprojection calculation was applied to Equation (A3) to obtain f(r), as: $\begin{matrix} f (r) = \int_{0}^{π} q (r sin θ) d θ \\ = 2 \int_{0}^{\infty} zP (z) \int_{0}^{π} cos (2 π zr sin θ) d θ dz \\ = 2 π \int_{0}^{\infty} zP (z) J_{0} (2 π rz) dz . \end{matrix}$ (A4)

The Bessel’s integral formula for J₀ (2πrz) [32] was used in the conversion. The substitution of Equation (A2) into Equation (A4) then yields: $f (r) = 2 π \sum_{n = 1}^{N} C_{n} {(\frac{2}{\sqrt{π}})}^{n} R^{\frac{n + 1}{2}} Γ (\frac{n}{2} + 1) \int_{0}^{\infty} z^{\frac{1 - n}{2}} J_{\frac{n + 1}{2}} (2 π Rz) J_{0} (2 π rz) dz .$ (A5)

In the case that 0 ≤ r < R, the integral formula (11.4.41) of [33] can be applied to Equation (A5) to obtain Equation (3). However, at r = R, different formulae were used as follows. According to the formula (6.574) of [31], all the terms for n≥2 on the right-hand side of Equation (A5) is zero. Only the term for n = 1 survives in Equation (A5), and the formula (11.4.42) of [33] was applied to the term to obtain Equation (4).

According to Equation (3) for 0≤r < R, the limiting behavior at the sample rim is that f (r) ⟶ C₁ as r ⟶R. Thus, it should be noted that f(r) ≈C₁ in the immediate vicinity of the sample rim, and is therefore independent of the N value (see Fig. 3) and the sample size (R). Because f(R) =C₁/2 (see Equation (4) and Fig. 2d), there is a gap in the f (r) values between f (r) as r ⟶R and f (R), which often occurs in the case of a Fourier transform involving a discontinuity.

The following equations may be helpful to the practical numerical calculation. We introduce the coefficient G_n that is independent of the sample chemistry and the X-ray source spectrum as defined by: $G_{n} = \frac{2}{\sqrt{π}} \frac{Γ (\frac{n}{2} + 1)}{Γ (\frac{n + 1}{2})} for n = 1, 2, 3, \dots, N,$ (A6)

to rewrite Equation (3) as: $f (r) = \sum_{n = 1}^{N} {2^{n - 1} G_{n} C_{n} (R^{2} - r^{2})}^{\frac{n - 1}{2}} for 0 \leq r < R .$ (A7)

The ratio of the gamma functions in Equation (A6) can be readily simplified to obtain: $G_{n} = \frac{2}{π} \frac{n}{G_{n - 1}},$

and $G_{0} = \frac{2}{π} for n = 1, 2, 3, \dots, N .$ (A8)

The use of the compact recurrence relation, Equation (A8), may be helpful to the effective computer programing for the calculation of f(r).

Appendix 2: The relationship between the non-linear projection and the nth moment of the LAC spectrum

The derivation of Equation (10) from Equations (8) and (9) is described herein. First, we derive recurrence relations on the nth moment by differentiation with respect to s. Second, we combine the recurrence relations with Equation (2) to relate μ_n(0) to C_n. Starting from Equations (7) to (9), we can derive the following set of equations: $M_{0} (s) = I (s),$ $\frac{{dM}_{n} (s)}{ds} = - M_{n + 1} (s) for n = 0, 1, 2, \dots, N - 1,$ $\frac{{dM}_{0} (s)}{ds} = \frac{dI}{ds} = - M_{1} (s),$ $μ_{n} (0) = \int_{0}^{\infty} i_{0} (E) μ (E)^{n} dE for n = 0, 1, 2, \dots, N,$ $μ_{0} (s) = 1,$ and $\begin{matrix} \frac{d μ_{n}}{ds} = \frac{1}{I^{2}} (\frac{{dM}_{n}}{ds} I - M_{n} \frac{dI}{ds}) = \frac{- M_{n + 1} I + M_{n} M_{1}}{I^{2}} = μ_{1} μ_{n} - μ_{n + 1} \\ for n = 0, 1, 2, \dots, N - 1 . \end{matrix}$ (A9)

Using the expressions in Equation (A9) with respect to μ_n, we obtain: $μ_{n + 1} = \sum_{m = 0}^{n} {(- 1)}^{m} (\begin{matrix} n \\ m \end{matrix}) μ_{n - m} \frac{d^{m} μ_{1}}{d s^{m}} for n = 1, 2, 3, \dots, N - 1,$ (A10)

where $\frac{d^{0} μ_{1}}{{ds}^{0}} = μ_{1},$ (A11)

and $(\begin{matrix} n \\ m \end{matrix}) = \frac{n!}{m! (n - m)!} .$ (A12)

Equation (A10) can be rewritten as: $\frac{d^{n} μ_{1}}{{ds}^{n}} = {(- 1)}^{n} μ_{n + 1} - \sum_{m = 0}^{n - 1} {(- 1)}^{n - m} (\begin{matrix} n \\ m \end{matrix}) μ_{n - m} \frac{d^{m} μ_{1}}{{ds}^{m}} for n = 1, 2, 3, \dots, N - 1 .$ (A13)

The moment described can be related to the power-series projection, h(s). Based on Equations (1), (9), and (A9), the first derivative of h(s) is expressed as: $\frac{dh}{ds} = \frac{dI}{ds} \frac{d}{dI} (- ln I) = \frac{M_{1}}{I} = μ_{1} .$ (A14)

The substitution of Equation (2) into Equation (A14) then yields: ${\frac{d^{n} h}{{ds}^{n}} |}_{s = 0} = {\frac{d^{n - 1} μ_{1}}{{ds}^{n - 1}} |}_{s = 0} = n! C_{n} for n = 1, 2, 3, \dots, N .$ (A15)

As a final step, Equation (A15) is substituted into Equation (A13) for s = 0 to relate μ_n(0) to C_n using the following recurrence relation:

$\begin{matrix} C_{n + 1} = - [\frac{{(- 1)}^{n + 1} μ_{n + 1} (0)}{(n + 1)!} + \frac{1}{(n + 1)!} \sum_{m = 0}^{n - 1} \frac{{(- 1)}^{n - m} n! (m + 1)! μ_{n - m} (0)}{m! (n - m)!} C_{m + 1}] \\ = - [\frac{{(- 1)}^{n + 1} μ_{n + 1} (0)}{(n + 1)!} + \sum_{m = 0}^{n - 1} \frac{{(- 1)}^{n - m} μ_{n - m} (0)}{(n - m)!} \frac{m + 1}{n + 1} C_{m + 1}] \\ for n = 0, 1, 2, 3, \dots, N - 1 . \end{matrix}$ (A16)

Equation (A16) can be simplified to Equation (10) using ν_n as defined by Equation (11). It should be noted that if n = 0, according to the conventional rule, the sum operation of the second term of the right-hand side of Equations (10) and (A16) is not performed and only the first term is calculated.

It may be helpful to explicitly show C_n for the first three values. These are: $C_{1} = μ_{1} (0),$ $C_{2} = \frac{μ_{1} (0)^{2} - μ_{2} (0)}{2},$ and $C_{3} = \frac{2 μ_{1} (0)^{3} - 3 μ_{1} (0) μ_{2} (0) + μ_{3} (0)}{6} .$ (A17)

The following simple relationship between the normalized X-ray intensity after sample penetration (i.e., I(s)/I(0)) and ν_n should be noted: $\frac{I (s)}{I (0)} = \int_{0}^{\infty} i_{0} (E) exp (- μ (E) s) dE = \int_{0}^{\infty} i_{0} (E) \sum_{n = 0}^{\infty} \frac{{[- μ (E) s]}^{n}}{n!} dE = \sum_{n = 0}^{\infty} ν_{n} s^{n} .$ (A18)

The power-series (i.e., Taylor-series) expansion of the natural exponential function was employed in the derivation of Equation (A18).

Footnotes

Acknowledgments

Comments by two anonymous reviewers were helpful. Preliminary CT experiments with cylindrical samples were performed at the GSJ-Lab AIST, the Japan Synchrotron Radiation Research Institute (proposal no. 2001B0501-NOD-np), and the Center for Advanced Marine Core Research (CMCR), Kochi University (13B034) with the support of JAMSTEC. The present study was supported in part by a KAKENHI Grant-in-Aid (no. 23241012) from the Japan Society for the Promotion of Science (JSPS).

References

Hammersberg

, Mångård, M., Correction for beam hardening artefacts in computerised tomography, J X-Ray Sci Technol 8 (1998), 75–93.

Kachelrieß

, Sourbelle

and Kalender

W.A.

, Empirical cupping correction: A first-order raw data precorrection for cone-beam computed tomography, Med Phys 33 (2006), 1269–1274.

Remeysen

, Swennen

, Beam hardening artifact reduction in microfocus computed tomography for improved quantitative coal characterization, Int J Coal Geol 67 (2006), 101–111.

Nakano

, Nakashima

, Nakamura

, Ikeda

, Observation and analysis of internal structure of rock using X-ray CT, J Geol Soc Jpn 106 (2000) 363–378 (in Japanese with English abstract).

Meganck

J.A.

, Kozloff

K.M.

, Thornton

M.M.

, Broski

S.M.

, Goldstein

S.A.

, Beam hardening artifacts in micro-computed tomography scanning can be reduced by X-ray beam filtration and the resulting images can be used to accurately measure BMD, Bone 45 (2009), 1104–1116.

Pietsch

, Jost

, Frenzel

, Raschke

, Walter

, Schirmer

, Hútter

, and Sieber

M.A.

, Efficacy and safety of lanthanoids as X-ray contrast agents, Eur J Radiol 80 (2011), 349–356.

Van de Casteele

, Van Dyck

, Sijbers

and Raman

, A model-based correction method for beam hardening artefacts in X-ray microtomography, J X-ray Sci Technol 12 (2004), 43–57.

Meng

, Zhang

, Wang

, Virtual experimentation of beam hardening effect in X-ray CT measurement of multiphase flow, Powder Technol 194 (2009), 153–157.

Nakashima

, The use of sodium polytungstate as an X-ray contrast agent to reduce the beam hardening artifact in hydrogeological laboratory experiments, J Hydrol Hydromech 61 (2013), 347–352.

10.

Nakashima

, Nakano

, Optimizing contrast agents with respect to reducing beam hardening in nonmedical X-ray computed tomography experiments, J X-Ray Sci Technol 22 (2014), 91–103.

11.

Park

H.S.

, Chung

Y.E.

, Seo

J.K.

, Computed tomographic beam-hardening artefacts: Mathematical characterization and analysis, Phil Trans R Soc A 373 (2015), 20140388.

12.

Arai

, Misawa

, Arai

, Shinozaki

, Sakamoto

, Yajima

, Nozaki

, Tajima

, Sato

, Hinoshita

, Accuracy analysis of intrahepatic fat density measurements using dual-energy computed tomography: Validation using a test phantom, J X-Ray Sci Technol 25 (2017), 403–415.

13.

Amato

, Asero

, Leotta

, Auditore

, Salamone

, Mannino

, Privitera

, Gueli

, Influence of the X-ray beam quality on the dose increment in CT with iodinated contrast medium, J X-Ray Sci Technol 24 (2016), 267–278.

14.

Nakashima

, Diffusivity measurement of heavy ions in Wyoming montmorillonite gels by X-ray computed tomography, J Contam Hydrol 61 (2003), 147–156.

15.

Jovanović

, Khan

, Enzmann

, Kersten

, Simultaneous segmentation and beam-hardening correction in computed microtomography of rock cores, Comp Geosci 56 (2013), 142–150.

16.

Herman

G.T.

Image Reconstruction from Projections, Academic Press, New York, 1980.

17.

Deans

S.R.

, Radon and Abel transforms (Chapter 8), Poularikas

A. D.

(ed.), The Transform and Applications Handbook (second edition), CRC Press, Boca Raton, 2000.

18.

Duan

, Wang

, Yu

, Leng

, McCollough

C.H.

, CT scanner X-ray estimation from transmission measurements, Med Phys 38 (2011), 993–997.

19.

Berger

M.J.

, Hubbell

J.H.

, Seltzer

S.M.

, Chang

, Coursey

J.S.

, Sukumar

, Zucker

D.S.

, Olsen

, XCOM: Photon Cross Sections Database (National Institute of Standards and Technology), http://www.nist.gov/pml/data/xcom/index.cfm/.

20.

Nakashima

, Nakano

, Nondestructive quantitative analysis of a heavy element in solution or suspension by single-shot computed tomography with a polychromatic X-ray source, Anal Sci 28 (2012), 1133–1138.

21.

Sidky

E.Y.

, Yu

, Pan

, Zou

, Vannier

, A robust method of x-ray source estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data, J Appl Phys 97 (2005), 124701.

22.

Zhao

, Niu

, Schafer

, Royalty

, An indirect transmission measurement-based estimation method for computed tomography, Phys Med Biol 60 (2015), 339–357.

23.

Lifton

J.J.

, Malcolm

A.A.

, McBride

J.W.

, A simulation-based study on the influence of beam hardening in X-ray computed tomography for dimensional metrology, J X-Ray Sci Technol 23 (2015), 65–82.

24.

Lifton

J.J.

, Multi-material linearization beam hardening correction for computed tomography, J X-ray Sci Technol 25 (2017), 629–640.

25.

Shi

, Yang

, Luo

, Reduce beam hardening artifacts of polychromatic X-ray computed tomography by an iterative approximation approach, J X-Ray Sci Technol 25 (2017), 417–428.

26.

Shah

J.P.

, Mann

S.D.

, Tornai

M.P.

, Characterization of X-ray scattering for various phantoms and clinical breast geometries using breast CT on a dedicated hybrid system, J X-Ray Sci Technol 25 (2017), 373–389.

27.

Tobin

K.W.

, Bingham

P.R.

, Gregor

, Mathematics of neutron imaging. In Neutron Imaging and Applications Anderson

I. S.

, McGreevy

R. L.

, and Billheux

H. Z.

(eds.) pp. 109–127, Springer, Boston, 2009.

28.

Zaidi

, Scheurer

A.H.

, Morel

, An object-oriented Monte Carlo simulator for 3D cylindrical positron tomographs, Comp Methods Prog Biomed 58 (1999), 133–145.

29.

Tanaka

H.K.

, Uchida

, Tanaka

, Shinohara

, Taira

, Development of a portable asbly-type cosmic-ray muon module for measuring the density structure of a column of magma, Earth, Planets Space 62 (2010), 119–129.

30.

De Rújula

, Glashow

S.L.

, Wilson

R.R.

and Charpak

, Neutrino exploration of the Earth, Phys Rep 99 (1983), 341–396.

31.

Gradshteyn

I.S.

, Ryzhik

I.M.

Table of Integrals, Series, and Products (sixth edition), Academic press, New York, 2000.

32.

Bowman

Introduction to Bessel Functions, Longmans, Green and Co., London, 1938.

33.

Abramowitz

, Stegun

I.A.

Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1972.

34.

Ketcham

R.A.

, Hanna

R.D.

, Beam hardening correction for X-ray computed tomography of heterogeneous natural materials, Comp Geosci 67 (2014), 49–61.

35.

Mail

, Moseley

D.J.

, Siewerdsen

J.H.

, Jaffray

D.A.

, The influence of bowtie filtration on cone-beam CT image quality, Med Phys 36 (2009), 22–32.

36.

Baba

, Konno

, Ueda

, Ikeda

, Comon of flat-panel detector and image-intensifier detector for cone-beam CT, Comput Med Imag Graph 26 (2002), 153–158.

37.

Goertzen

A.L.

, Beekman

F.J.

, Cherry

S.R.

, Effect of phantom voxelization in CT simulations, Med Phys 29 (2002), 492–498.