Non-convex optimization based optimal bone correction for various beam-hardening artifacts in CT imaging

Abstract

Tube of X-ray computed tomography (CT) system emitting a polychromatic spectrum of photons leads to beam hardening artifacts such as cupping and streaks, while the metal implants in the imaged object results in metal artifacts in the reconstructed images. The simultaneous emergence of various beam-hardening artifacts degrades the diagnostic accuracy of CT images in clinics. Thus, it should be deeply investigated for suppressing such artifacts. In this study, data consistency condition is exploited to construct an objective function. Non-convex optimization algorithm is employed to solve the optimal scaling factors. Finally, an optimal bone correction is acquired to simultaneously correct for cupping, streaks and metal artifacts. Experimental result acquired by a realistic computer simulation demonstrates that the proposed method can adaptively determine the optimal scaling factors, and then correct for various beam-hardening artifacts in the reconstructed CT images. Especially, as compared to the nonlinear least squares before variable substitution, the running time of the new CT image reconstruction algorithm decreases 82.36% and residual error reduces 55.95%. As compared to the nonlinear least squares after variable substitution, the running time of the new algorithm decreases 67.54% with the same residual error.

Keywords

X-ray CT;beam hardening;consistency condition;non-convex optimization

1 Introduction

Tube of X-ray CT imaging system emits photons with polychromatic spectrum, resulting in beam hardening artifacts in the reconstructed image, such as cupping and streaks artifacts [1 –3]. At the same time, there may also be metal implants in the imaged object, resulting in metal artifacts in the reconstructed image. For a long time, metal artifacts correction is a difficult problem in X-ray CT imaging of modern medical diagnosis. The simultaneous emergence of various types of artifacts certainly will degrade the accuracy of clinical diagnosis, so it is necessary to conduct in-depth research to suppress them.

Traditional water correction can only correct for cupping artifacts caused by water like substances [4, 5]. Traditional bone correction can correct for the cupping artifact caused by water like substances and the streaks artifacts caused by bone like substances at the same time, but it is usually necessary to use cumbersome physical experiments to determine the scale factor, λ₀, in advance, which is subsequently used to convert the reconstruction value of bone tissue into bulk density [6, 7]. We have ever proposed an optimal bone correction method based on the data consistency condition, which can determine the optimal scale factor [8]. The optimal bone correction method can not only avoid large number of physical experiments, but also adapt to the changes of imaging equipment and imaging objects.

Based on the mathematical and physical model of X-ray imaging [9], this paper constructed the corresponding objective function by using data consistency condition for the case that the imaged objects include water-like, bone-like and metal-like substances at the same time. Essentially, the objective function constructed on the novel mathematical and physical model is not globally convex, therefore, the alternating optimization algorithm [10] and non-linear least squares cannot guarantee convergence unless the initial value of variables is very close to the global optimal point. Therefore, non-convex optimization algorithm [11] was used to adaptively determine the optimal scale factors of soft tissue, bone and metal, respectively. An optimal bone correction method was then obtained to correct for the cupping, streaks and metal artifacts simultaneously. Finally, the performance of the correction method is verified by computer simulation.

The innovations of the proposed method are at least three aspects: (1) the beam-hardening effects of the metal-like substance is explicitly considered in the mathematical and physical model of X-ray imaging; (2) the superiority of non-convex optimization algorithm and (3) the positive influence of variable substitution on the performance of the optimal bone correction for various beam-hardening artifacts in x-ray CT imaging are validated in this study.

2 Methodology

2.1 Mathematical model of X-ray imaging

The data acquisition process of projections, generated by X-ray source with polychromatic spectrum in medical diagnostic CT system, can be mathematically modeled as:

${\hat{g}}_{i}^{f} = log (\frac{I_{0, i}^{poly}}{I_{i}^{poly}}) = log (\frac{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) d ɛ}{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) exp (- \int_{L_{i}} μ (l, ɛ) ρ (l) dl) d ɛ})$ (1) where the functions S_i (ɛ) and Q_i (ɛ) with independent variable ɛ represent the emission spectrum and the absorption spectrum distributions of the source and detector of i^th X-ray, respectively, and the data consistency condition described in appendix is not satisfied between the projection data set ${{\hat{g}}_{i}^{f}}$ and image function f^mono (x, y) = μ (x, y, ɛ_mono) ρ (x, y).

Analytical CT reconstruction usually assumes that the photons emitted from the X-ray source used in projection data acquisition are of monochromatic energy. Assuming that the monochromatic energy of the photons emitted by X-ray source is ɛ_mono, the corresponding projection data acquisition process can be mathematically modeled as

$g_{i}^{f} = log (\frac{I_{0, i}^{mono}}{I_{i}^{mono}}) = log (\frac{S_{i} (ɛ_{mono}) Q_{i} (ɛ_{mono})}{S_{i} (ɛ_{mono}) Q_{i} (ɛ_{mono}) exp (- \int_{L_{i}} μ (l, ɛ_{mono}) ρ (l) dl)}) = \int_{L_{i}} μ (l, ɛ_{mono}) ρ (l) dl$ (2) where the data consistency conditions described in appendix are satisfied between projection data set ${g_{i}^{f}}$ and image function f^mono (x, y) = μ (x, y, ɛ_mono) ρ (x, y).

The beam hardening artifacts in the reconstructed image can be thought to be caused by the mismatch between the assumption of the monochromatic spectrum in the analytical CT reconstruction and the fact of the polychromatic spectrum of the projected data set.

2.2 Water correction

Water correction is a traditional beam hardening correction method suitable for single material [4, 5], which can be divided into two steps: calibration and correction.

In the calibration step, water-like substances (such as plexiglass plates of different thicknesses) are taken as the imaged objects to obtain the projection data of polychromatic spectrum and determine the coefficients of the correction polynomial as follows:

$g = \sum_{n = 0}^{N} a_{n} {\hat{g}}^{n}$ (3)

In the correction step, the projection value of monochromatic spectrum corresponding to the actual projection data ${\hat{g}}_{i}^{f}$ was calculated according to the above correction polynomial

$g_{i}^{f} = P^{w} [{\hat{g}}_{i}^{f}] ≜ \sum_{n = 0}^{N} a_{n} {({\hat{g}}_{i}^{f})}^{n}$ (4)

for subsequent analytical CT reconstruction. The effectiveness of water correction depends on whether the imaged object is water-like.

2.3 Bone correction

When bone correction [6, 7] is performed, the projection is firstly water corrected [4, 5], and then analytic CT pre-reconstruction is performed to obtain water corrected image f^w (x, y). The soft tissue and bone tissue images were obtained by thresholding the pre-reconstructed image,

$f^{s} (x, y) ≜ {\begin{matrix} f^{w} (x, y), f^{w} (x, y) < T \\ 0, f^{w} (x, y) ⩾ T \end{matrix}$ (5)

$f^{b} (x, y) ≜ {\begin{matrix} 0, f^{w} (x, y) < T \\ f^{w} (x, y), f^{w} (x, y) ⩾ T \end{matrix}$ (6) where the parameter T indicates the manual threshold. Then, the area densities of soft tissue and bone tissue were obtained by re-projection:

$ρ_{i}^{a, s} = \frac{f_{i}^{s}}{μ_{s} (ɛ_{mono})} = \frac{\iint_{(x, y) \in L_{i}} f^{s} (x, y) dxdy}{μ_{w} (ɛ_{mono})}$ (7)

$ρ_{i}^{a, b} = \frac{f_{i}^{b}}{μ_{b} (ɛ_{mono})} = \frac{\iint_{(x, y) \in L_{i}} f^{b} (x, y) dxdy}{μ_{w} (ɛ_{mono}) λ_{0}}$ (8) where μ_w, μ_s and μ_b are the mass attenuation coefficients of water, soft tissue and bone respectively, μ_s ≈ μ_w. The scale factor λ₀ is approximately equivalent to μ_b (ɛ_mono)/ μ_w (ɛ_mono), which is used to convert the bone tissue reconstruction value into the bone tissue area density. In the physical sense, there is 1 < λ₀.

Finally, the bone correction is formulized as follows,

$\begin{array}{l} g_{i}^{f} = P^{w} [{\hat{g}}_{i}^{f}] + f_{i} - P^{w} (\log (\frac{\int_{ε_{\min}}^{ε_{\max}} S_{i} (ε) Q_{i} (ε) d ε n}{\int_{ε_{\min}}^{ε_{\max}} S_{i} (ε) Q_{i} (ε) d ε \exp (- {\dot{μ}}_{s} f_{i}^{s} - {\dot{μ}}_{b} f_{i}^{b} / λ_{0}) d ε})) \\ {\dot{μ}}_{s} (ε) = μ_{s} (ε) / μ_{w} (ε_{m o n o}), {\dot{μ}}_{b} (ε) = μ_{b} (ε) / μ_{w} (ε_{m o n o}) \end{array}$ (9) where $f_{i} = f_{i}^{s} + f_{i}^{b}$ , and $P^{w} [{\hat{g}}_{i}^{f}] + f_{i}$ doesn’t change with λ₀.

2.4 Optimal bone correction

When metal exists in the imaged object, the pre-reconstructed image after water correction is segmented by thresholding to obtain soft tissue, bone tissue and metal images, respectively $f^{s} (x, y) ≜ {\begin{matrix} f^{w} (x, y), f^{w} (x, y) < T_{1} \\ 0, f^{w} (x, y) ⩾ T_{1} \end{matrix}$ (10) ∥ $f^{b} (x, y) ≜ {\begin{matrix} 0, f^{w} (x, y) < T_{1} \\ f^{w} (x, y), T_{1} ⩽ f^{w} (x, y) < T_{2} \\ 0, T_{2} ⩽ f^{w} (x, y) \end{matrix}$ (11) ∥ $f^{mt} (x, y) ≜ {\begin{matrix} 0, f^{w} (x, y) < T_{2} \\ f^{w} (x, y), T_{2} ⩽ f^{w} (x, y) \end{matrix}$ (12) where the parameters T₁ < T₂ are the manual threshold values. Then, the area densities of soft tissue, bone tissue and metal implant were respectively obtained by re-projecting their respective images as follows, $ρ_{i}^{a, s} = \frac{f_{i}^{s}}{μ_{s} (ɛ_{mono})} = \frac{\iint_{(x, y) \in L_{i}} f^{s} (x, y) dxdy}{μ_{w} (ɛ_{mono}) λ_{- 1}}$ (13) ∥ $ρ_{i}^{a, b} = \frac{f_{i}^{b}}{μ_{b} (ɛ_{mono})} = \frac{\iint_{(x, y) \in L_{i}} f^{b} (x, y) dxdy}{μ_{w} (ɛ_{mono}) λ_{0}}$ (14) ∥ $ρ_{i}^{a, mt} = \frac{f_{i}^{mt}}{μ_{mt} (ɛ_{mono})} = \frac{\iint_{(x, y) \in L_{i}} f^{mt} (x, y) dxdy}{μ_{w} (ɛ_{mono}) λ_{1}}$ (15) where μ_mt is the metal mass attenuation coefficient. The scale factor λ_-1 is approximately equivalent to μ_s (ɛ_mono)/ μ_w (ɛ_mono) and is used to convert the soft tissue reconstruction value into the soft tissue area density. The scale factor λ₁ is approximately equivalent to μ_mt (ɛ_mono)/ μ_w (ɛ_mono) and is used to convert the metal reconstruction value into the metal area density. In the physical sense, there is 1 < λ_-1 < λ₀ < λ₁.When metal exists in the scanned object, the bone correction formula (9) is changed to: $\begin{matrix} g_{i}^{f} = P^{w} [{\hat{g}}_{i}^{f}] + f_{i} - P^{w} (log (\frac{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) d ɛ}{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) exp (- {\dot{μ}}_{s} f_{i}^{s} / λ_{- 1} - {\dot{μ}}_{b} f_{i}^{b} / λ_{0} - {\dot{μ}}_{mt} f_{i}^{mt} / λ_{1}) d ɛ})) \\ {\dot{μ}}_{mt} (ɛ) ≜ μ_{mt} (ɛ) / μ_{w} (ɛ_{mono}) \end{matrix}$ (16) where $f_{i} = f_{i}^{s} + f_{i}^{b} + f_{i}^{mt}$ , and $P^{w} [{\hat{g}}_{i}^{f}] + f_{i}$ doesn’t change with λ_-1, λ₀ and λ₁.

2.4.1 Nonlinear least squares

When metal exists in the scanned object, the optimal scale factor of polychromatic spectrum in X-ray CT imaging can be determined by minimizing the following objective function by nonlinear least squares: $\begin{matrix} Φ_{1} (g^{f}, λ_{- 1}, λ_{0}, λ_{1}) = \sum_{φ} {\int_{- \infty}^{+ \infty} {g^{f} (t, φ - atan (t / D)) - \\ f (t, φ - atan (t / D)) - P^{w} ({\hat{g}}^{f} (t, φ - atan (t / D))) + \\ P^{w} [log (\frac{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) d ɛ}{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) exp (- {\dot{μ}}_{s} f_{i}^{s} / λ_{- 1} - {\dot{μ}}_{b} f_{i}^{b} / λ_{0} - {\dot{μ}}_{mt} f_{i}^{mt} / λ_{1}) d ɛ})]} {\frac{D^{3}}{{(D^{2} + t^{2})}^{1.5}} dt}}^{2} \end{matrix}$ (17) However, the unknowns g^f in Φ₁ is a function of (t, φ), and there are too many unknown variables, so it is difficult to realize the minimization process.Because the projection f = f^s + f^b + f^mt naturally satisfies the data consistency condition (see Appendix), and g^f can reduce the number of unknowns by using data consistency condition, Φ₁ can be further transformed into: $\begin{matrix} Φ_{2} (m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) = \sum_{φ} {(m_{0, 0} - 1) m_{0, 0}^{f} + \\ \int_{- \infty}^{+ \infty} {P^{w} [log (\frac{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) d ɛ}{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) exp (- {\dot{μ}}_{s} f_{i}^{s} / λ_{- 1} - {\dot{μ}}_{b} f_{i}^{b} / λ_{0} - {\dot{μ}}_{mt} f_{i}^{mt} / λ_{1}) d ɛ})] - \\ P^{w} [{\hat{g}}^{f} (t, φ - atan (t / D))]} {\frac{D^{3}}{{(D^{2} + t^{2})}^{1.5}} dt}}^{2} \end{matrix}$ (18) where $m_{0, 0}^{f}$ can be calculated in advance as follows: $m_{0, 0}^{f} = \int_{- \infty}^{+ \infty} f (t, φ - atan (t / D)) \frac{D^{3}}{{(D^{2} + t^{2})}^{1.5}} dt$ (19) and notice that $m_{0, 0}^{f}$ doesn’t depend on (t, φ).For the nonlinear least squares problem expressed in Equation (18), the alternating optimization algorithm [10] cannot guarantee convergence unless the initial value of variables is very close to the global optimal point. Therefore, this paper tries to solve this problem through non-convex optimization.

2.4.2 Non-convex optimization

Since the objective function Φ₂ (see Equation (18)) cannot guarantee the global convexity, the non-convex optimization algorithm [11] should be used to solve the problem on the premise of reasonably selecting the initial value of variables $(m_{0, 0}^{(0)}, λ_{- 1}^{(0)}, λ_{0}^{(0)}, λ_{1}^{(0)})$ , and that Equation (16) is replaced by the simplified nonlinear model as follows, $\begin{matrix} G (φ) ≜ \int_{- \infty}^{+ \infty} P^{w} [{\hat{g}}^{f} (t, φ - atan (t / D))] \frac{D^{3}}{{(D^{2} + t^{2})}^{1.5}} dt \\ \approx G^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) ≜ (m_{0, 0} - 1) m_{0, 0}^{f} + \\ \int_{- \infty}^{+ \infty} P^{w} [log (\frac{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) d ɛ}{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) exp (- {\dot{μ}}_{s} f_{i}^{s} / λ_{- 1} - {\dot{μ}}_{b} f_{i}^{b} / λ_{0} - {\dot{μ}}_{mt} f_{i}^{mt} / λ_{1}) d ɛ})] \frac{D^{3}}{{(D^{2} + t^{2})}^{1.5}} dt \end{matrix}$ (20) Then, the procedure of the non-convex optimization is listed as the pseudo code below:(1) Selecting the initial value of variables: ${\tilde{m}}_{0, 0} = m_{0, 0}^{(0)}, {\tilde{λ}}_{- 1} = λ_{- 1}^{(0)}, {\tilde{λ}}_{0} = λ_{0}^{(0)}, {\tilde{λ}}_{1} = λ_{1}^{(0)})$ ;(2) Taylor expansion of Equation (20) at $m_{0, 0} = {\tilde{m}}_{0, 0}, λ_{- 1} = {\tilde{λ}}_{- 1}, λ_{0} = {\tilde{λ}}_{0}, λ_{1} = {\tilde{λ}}_{1})$ : $\begin{matrix} G (φ) \approx G^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) = G^{NL} (φ, {\tilde{m}}_{0, 0}, {\tilde{λ}}_{- 1}, {\tilde{λ}}_{0}, {\tilde{λ}}_{1}) + m_{0, 0}^{f} (m_{0, 0} - {\tilde{m}}_{0, 0}) \\ + K_{λ_{- 1}} (φ) (λ_{- 1} - {\tilde{λ}}_{- 1}) + K_{λ_{0}} (φ) (λ_{0} - {\tilde{λ}}_{0}) + K_{λ_{1}} (φ) (λ_{1} - {\tilde{λ}}_{1}) + vartriangle G^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}), \end{matrix}$ (21) where $m_{0, 0}^{f}$ , K_{λ_-1} (φ), K_λ₀ (φ) and K_λ₁ (φ) are the partial derivatives of G^NL (φ, m_0,0, λ_-1, λ₀, λ₁) with respect to m_0,0, $λ_{- 1}$ , $λ_{0}$ and $λ_{1}$ at ${\tilde{m}}_{0, 0}$ , ${\tilde{λ}}_{- 1}$ , ${\tilde{λ}}_{0}$ and ${\tilde{λ}}_{0}$ , respectively, and vartriangleG^NL is an infinitesimal of higher order. ${\tilde{m}}_{0, 0}$ , ${\tilde{λ}}_{- 1}$ , ${\tilde{λ}}_{0}$ and ${\tilde{λ}}_{1}$ are the results at the previous iteration step;(2.1) Similar to the non-convex optimization algorithm [11], the following nonlinear least squares problem is solved iteratively: $G (φ) - vartriangle {\tilde{G}}^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) = [\begin{matrix} m_{0, 0}^{f} & K_{λ_{- 1}} (φ) & K_{λ_{0}} (φ) & K_{λ_{1}} (φ) \end{matrix}] [\begin{matrix} m_{0, 0} & λ_{- 1} & λ_{0} & λ_{1} \end{matrix}]$ (22) where $\begin{matrix} vartriangle {\tilde{G}}^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) ≜ G^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) \\ - [\begin{matrix} m_{0, 0}^{f} & K_{λ_{- 1}} (φ) & K_{λ_{0}} (φ) & K_{λ_{1}} (φ) \end{matrix}] [\begin{matrix} m_{0, 0} & λ_{- 1} & λ_{0} & λ_{1} \end{matrix}] \end{matrix}$ (23) Consequently, ${\tilde{m}}_{0, 0}$ , ${\tilde{λ}}_{- 1}$ , ${\tilde{λ}}_{0}$ and ${\tilde{λ}}_{1}$ are updated;(3) After the optimal scale factors $(λ_{- 1}^{*}, λ_{0}^{*}, λ_{1}^{*})$ are determined, the artifacts correction is carried out according to the bone correction formula (16).

2.4.3 Variable substitution

The simplified model (Equation (20)) is nonlinear, and the objective function Φ₂ cannot guarantee global convexity. In order to simplify the solving process of the non-convex optimization algorithm [11], the scale factors $(λ_{- 1}, λ_{0}, λ_{1})$ are replaced by $(λ_{- 1}^{- 1}, λ_{0}^{- 1}, λ_{1}^{- 1})$ , thus obtaining the following nonlinear simplified model: $\begin{matrix} G^{NL} (φ, m_{0, 0}, λ_{- 1}, λ_{0}, λ_{1}) ≜ (m_{0, 0} - 1) m_{0, 0}^{f} + \\ \int_{- \infty}^{+ \infty} P^{w} [log (\frac{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) d ɛ}{\int_{ɛ_{min}}^{ɛ_{max}} S_{i} (ɛ) Q_{i} (ɛ) exp (- {\dot{μ}}_{s} f_{i}^{s} λ_{- 1} - {\dot{μ}}_{b} f_{i}^{b} λ_{0} - {\dot{μ}}_{mt} f_{i}^{mt} λ_{1}) d ɛ})] \frac{D^{3}}{{(D^{2} + t^{2})}^{1.5}} dt \end{matrix}$ (24) and after the optimal scale factors $(λ_{- 1}^{*}, λ_{0}^{*}, λ_{1}^{*})$ are determined, the inverse substitution was $((λ_{- 1}^{*})^{- 1}, (λ_{0}^{*})^{- 1}, (λ_{1}^{*})^{- 1})$ , and then the artifacts correction is carried out according to the bone correction formula (16).

Note that this variable substitution is empirically beneficial to improve the linearity of the simplified model formula (20) and the global convexity of the objective function Φ₂.

3 Computer simulation verification

Our previous approach to analytical projection calculation [9] of FORBILD head phantom [12] at z = 0mm section is used for the computer simulation verification. The computational resources used to conduct the simulation are Windows7 OS and Intel^® Core™ i3-2310M CPU @2.10 GHz.

3.1 Mechanical parameters

In the simulation experiment, the radius of CT circular scanning trajectory is set as 50cm. Fig. 1 shows the natural coordinate system, the local coordinate system of the two-dimension X-ray detector and the scanning track of the X-ray source point. The center of the phantom coincides with the rotating center O, and the radius of the phantom is 128mm. The X-ray detector contains pixels, which constitute an 850×80 matrix, and each pixel size is 1mm×1mm. The distance between the point O_d and the detector center is 50cm. When the X-ray source point is rotated 360°, the projection data of 1080 view angles are evenly collected.

Fig. 1

Schematic of natural coordinate system, local coordinate system of X-ray 2D detector, and scanning trajectory of X-ray source point.

3.2 Energy spectrum distributions

The performance of the proposed method is independent on the types of X-ray tube and detector. X-ray tube emission spectrum S (ɛ) [13] and CsI material X-ray detector absorption spectrum Q (ɛ) [14] are used and shown in Fig. 2, in which the tube voltage used for polychromatic spectrum imaging is set to 120kVp (plus 35mm aluminum filter), and the monochromatic spectrum used for comparison is manually set at 66keV. Under this monochromatic spectrum, equivalent linear attenuation coefficient close to polychromatic spectrum can be obtained. It is convenient for follow-up experiment comparison.

Fig. 2

Distribution functions of emitting energy spectrum S (ɛ) as X-ray tube is set to 120 kVp and a filter of 35 mm Aluminum is added, and absorption energy spectrum Q (ɛ) of CsI X-ray detector.

3.3 Analytical CT reconstruction

Analytical CT reconstruction adopts the algorithm in paper [15], which is a typical filtering back projection (FBP) algorithm. The reconstructed image is a matrix of 512×512 pixels, with each pixel size of 0.5mm×0.5mm. Since the two-dimensional reconstruction plane locates in the middle plane of the cone beam, the reconstruction is theoretically accurate. Fig. 3 shows the reconstruction results of monochromatic spectrum and polychromatic spectrum, in which a small bubble in the right ear has been intentionally replaced by metal material iron to represent the beam-hardening and metal artifacts caused by a small volume metal implant.

3.4 Optimal bone correction

Fig. 3

Reconstructed slice images of FORBILD head phantom at z = 0mm section. The display window of the reconstructed images is with [–5.4 114.1] HU, while that of the difference images with [–54 54] HU.

The optimal bone correction adopts non-convex optimization [11] after variable substitution and is compared with nonlinear least squares. Since the proposed minimization problem cannot guarantee global convexity, in order to compare the performance of each algorithm, three different groups of the initial values of variables are selected: (i) $m_{0, 0}^{(0)} = 1, λ_{- 1}^{(0)} = 1, λ_{0}^{(0)} = 1, λ_{1}^{(0)} = 4.5$ , (ii) $m_{0, 0}^{(0)} = 1, λ_{- 1}^{(0)} = 10, λ_{0}^{(0)} = 10, λ_{1}^{(0)} = 10$ , (iii) $m_{0, 0}^{(0)} = 1, λ_{- 1}^{(0)} = 0.1, λ_{0}^{(0)} = 0.1, λ_{1}^{(0)} = 0.1$ . When non-convex optimization [11] is adopted, no matter what initial value is selected, there is always the optimal convergence result at the end of iteration: $m_{0, 0}^{*} = 1.1831$ , $λ_{- 1}^{*} = 1.2449$ , $λ_{0}^{*} = 1.5068$ , $λ_{1}^{*} = 5.0281$ .

Figs 4 -9 show the changing curves of the iterative process of non-convex optimization [11] under the condition that the first group of the initial values of variables is selected.

Fig. 4

Changing tendency of residual error Φ₂ during the optimal bone correction with non-convex optimization.

Fig. 5

Profile of Fig.4 shown using semiology.

Fig. 6

Changing tendency of λ_–1 during the optimal bone correction with non-convex optimization.

Fig. 7

Changing tendency of λ₀ during the optimal bone correction with non-convex optimization.

Fig. 8

Changing tendency of λ₁ during the optimal bone correction with non-convex optimization.

Fig. 9

Changing tendency of m₀₀ during the optimal bone correction with non-convex optimization.

3.5 Data analysis

Fig. 10 qualitatively investigates various correction effects, in which Equation (9) is adopted for bone correction without considering the effect of λ₁. As can be seen from the results in Fig. 10, the cupping artifacts are easily corrected by all methods. The soft tissue areas adjacent to the right ear of FORBILD head phantom were selected to examine the streak artifacts, which are caused by a combination of bone tissue and metal implant.

Fig.10

Result comparison between different optimization methods after variable substitution. Reconstructed slice images of FORBILD head phantom at z = 0 mm section are displayed with window [–5.4 114.1] HU, while the difference images are with [–54 54] HU

Water correction has little effect on the streak artifacts, and the streaks severity is almost the same as that of the polychromatic spectrum reconstruction. Bone correction (λ₀ = 1.0) showed overcorrection, the dark streaks caused by bone tissue became brighter, and the cupping artifacts was overcorrected to capping artifacts. Bone correction (λ₀ = 1.6) showed under-correction, and there were still dark streaks caused by bone tissue. Bone correction (λ₀ = 1.3) as it is closest to the optimal bone correction, the correction effect is very close. However, all bone corrections didn’t take into account the metal effects, and the metal artifacts caused by the small metal implant were not corrected for, presenting as residual dark streaks. As expected, the optimal bone correction (using Equation (16)) performed well with almost no metal artifact residues. As the comparison, the metal artifact correction method with a linear interpolation in the projection domain is also carefully performed [16]. Due to the large amount of high-frequency content in the right ear area of the FORBILD head phantom, the metal artifact correction performance degrades seriously.

Figs 11 12 are the profiles of various correction reconstruction and monochromatic and polychromatic spectrum reconstruction images. The studied line is at the line 385 of the FORBILD head phantom z = 0mm section, while shown in Fig. 12 are the profiles of Fig. 11 after zoomed in.

Fig. 11

Profiles of various corrected and monochromatic reconstructions along the 385th line in the images of FORBILD head phantom at z = 0mm section.

Fig. 12

Profiles of Fig. 11 after zoomed in.

In order to further quantitatively analyze and compare the correction effects, ${\tilde{σ}}_{b}$ and ${\tilde{σ}}_{s}$ are introduced to represent the degrees of inconsistencies in the bone and soft tissue regions, respectively,

${\tilde{σ}}_{b, s} = σ_{b, s} / mea n_{b, s} \times 100 %$ (25) where

$mea n_{b, s} = {\sum_{(x, y) \in R_{b, s}} f^{c} (x, y)} / N_{b, s}$ (26)

$σ_{b, s} = {(\sum_{(x, y) \in R_{b, s}} {(f^{c} (x, y) - mea n_{b, s})}^{2} / (N_{b, s} - 1))}^{0.5}$ (27) and f^c (x, y) is the reconstructed image. R_b and R_s represent bone and soft tissue regions, respectively. N_b and N_s represent the total numbers of pixels in bone and soft tissue regions, respectively. Bone and soft tissue regions were defined as the regions whose CT values are 50 Hu and 800Hu in the original definition of FORBILD head phantom, respectively. Meanwhile, ${\tilde{δ}}_{b}$ and ${\tilde{δ}}_{s}$ are introduced to represent the deviations from monochromatic result in the bone and soft tissue regions, respectively,

${\tilde{δ}}_{b, s} = δ_{b, s} / mea n_{b, s} \times 100 %$ (28) where

$δ_{b, s} = {(\sum_{(x, y) \in R_{b, s}} {(f^{c} (x, y) - f^{mono} (x, y))}^{2} / (N_{b, s} - 1))}^{0.5}$ (29)

The quantitative analysis results in Table 1 show that the optimal bone correction result does exist, and the non-convex optimal bone correction can balance the correction performance in the different regions well, and show an excellent correction effect on metal artifacts (See Fig. 10).

Table 1

Inconsistencies and deviations in soft tissue region, whose CT value is 50Hu, and bone region, whose CT value is 800Hu, in the original definition of FORBILD head phantom. Unit: %

	Inconsistency		Deviation
Method	Bone tissue	Soft tissue	Bone tissue	Soft tissue
Monochromatic	1.6251	2.0158	0	0
Polychromatic	2.0190	1.8914	1.3743	0.5661
Water correction	1.8349	1.8649	1.0576	0.4367
Bone correction (λ_0 =1.0)	2.6898	2.0432	6.1647	0.4731
Bone correction (λ_0 =1.3)	1.5790	1.8528	0.4670	0.2665
Bone correction (λ_0 =1.6)	1.0728	1.8021	1.0576	0.3601
Non-convex optimal bone correction	1.5162	1.7322	0.4060	0.3650
Metal artifact correction	2.7556	1.7918	3.0507	2.4614

In order to verify the necessity of non-convex optimization for optimal bone correction, the third group of the initial values of variables was selected. Non-convex optimization [11] before/after variable substitution is compared with nonlinear least squares, and the comparison results are shown in Fig. 13 and Table 2. It can be seen from Table 2, that the optimal bone correction with non-convex optimization after variable substitution has the minimum running time and residual error Φ₂, indicating that both variable substitution and non-convex optimization are very necessary. Especially, as compared to nonlinear least squares before variable substitution, the running time decreases 82.36% and residual error Φ₂ reduces 55.95%. As compared to nonlinear least squares after variable substitution, the running time decreases 67.54% with the same residual error Φ₂.

Fig. 13

Result comparison of the optimal bone corrections with nonlinear least squares and non-convex optimization before (the first row) and after (the last row) variable substitution, wherein $m_{0, 0}^{(0)} = 1, λ_{- 1}^{(0)} = 0.1, λ_{0}^{(0)} = 0.1, λ_{1}^{(0)} = 0.1$ . Reconstructed slice images of FORBILD head phantom at z = 0mm section are displayed with window [–22 141] HU.

Table 2

Performance comparison of the optimal bone corrections with nonlinear least squares and non-convex optimization before and after variable substitution, wherein $m_{0, 0}^{(0)} = 1, λ_{- 1}^{(0)} = 0.1, λ_{0}^{(0)} = 0.1, λ_{1}^{(0)} = 0.1$

Variable Substitution	Method	m ₀₀	λ_–1	λ₀	λ₁	Time/s	Φ ₂
Before	Nonlinear least squares	–1.4530	0.2641	0.5816	0.4072	8626.6	3039.7
	Non-convex optimization	1.1831	1.2449	1.5068	5.0281	1521.5	1949.1
After	Nonlinear least squares	1.1831	1.2449	1.5068	5.0281	2920.5	1949.1
	Non-convex optimization	1.1831	1.2449	1.5068	5.0281	947.93	1949.1

4 Discussions

Metal artifacts usually include multiple factors: beam hardening artifacts caused by a polychromatic spectrum, high-frequency streaks artifacts reconstructed from metal projection data with a large dynamic range that are tangent to the edge of metal implants, photon starvation, quantum noise and so on. In this paper, only soft tissue and bone induced beam hardening artifacts (cupping and dark streaks) are considered. Therefore, the beam hardening artifacts (dark streaks) caused by small metal implants can be effectively solved in the framework of optimal bone correction. Computer simulation experiments and analysis results also preliminarily verify the performance of the proposed method. High-frequency streaks artifacts correction is not within the scope of discussion in this paper, and filtering kernel parameters may be alleviated to a certain extent by adaptive adjustment reconstruction algorithm, which will be reserved for future research.

In this study, only metal implants with small volume are discussed. Metal implants with large volume will cause serious photon starvation and quantum noise, leading to serious degradation of projection data quality, and the optimal bone correction method proposed in this paper needs to be improved appropriately [17]. When it is difficult to obtain the energy spectrum distribution functions of X-ray source emission and detector absorption, the polychromatic spectrum effect can be approximated by polynomial. In this paper, only analytical simulation experiment with FORBILD head phantom is used to verify the proposed method, but as a result that the phantom has long been widely used, and the performance of CT reconstruction and correction algorithm has good performance, so the simulation results are still reliable to some extent, and we will try to collect actual data for the goal of verification.

From the experimental results, it can be observed that the optimal bone correction with non-convex optimization after variable substitution has the minimum running time and residual error, indicating that both variable substitution and non-convex optimization are very necessary. As for the underlying cause on the usefulness of variable substitution, it is still not completely clear. One possible explanation is that variable substitution can change the morphology of residual error Φ₂ to some extent.

In recent years, deep-learning (DL) technique has developed very quickly, and the effect of its applications to various engineering problems is significant. As DL is used in the correction for the various artifacts in X-ray CT imaging, the effect is obvious too [18, 19]. The relevant researches based on DL by combining with results reported in this paper will be our future works.

5 Conclusions

In this paper, a non-convex optimal bone correction method based on data consistency and suitable for both beam hardening artifacts and metal artifacts is proposed. This method can determine the optimal values of all scale factors adaptively. The results of computer numerical simulation and qualitative and quantitative analysis demonstrate the existence of the optimal scale factors of bone correction and the adaptive determination ability of the proposed non-convex optimal bone correction method to the scale factors.

Footnotes

Acknowledgments

The project is supported in part by the grants from Xi’an Key Laboratory of Advanced Controlling and Intelligent Processing, China (2019220714SYS022CG044), Science and Technology Department of Shaanxi Province (2020SF-377), China, and National Natural Science Foundation of China (N0. 62071378).

Appendix: Data Consistency Conditions

References

Zhang

, Li

, Zhang

. et al., Review of the methods for beam hardening correction in X-ray computed tomography, CT Theory and Applications 22(1) (2013), 195–204.

Chen

, Pan

A beam-hardening correction method based on poly-energetic statistics in X-ray CT, CT Theory and Applications 19(1) (2010), 21–27.

, Zhang

Estimation of energy spectra of CT system and beam-hardening correction algorithm, CT Theory and Applications 12(2) (2003), 10–15.

Herman

G.T.

Correction for beam hardening in computed tomography, Physics in Medicine and Biology 24(1) (1979), 81–106.

Tang

, Lin

Y.B.

, Sun

G.Y.

, Bao

X.D.

A metal artifact reduction scheme in CT by a Poisson fusion sinogram based postprocessing method, Journal of X-ray Science and Technology 29(2) (2021), 245–257.

Joseph

P.M.

, Spital

R.D.

A method for correcting bone induced artifacts in computed tomography scanners, Journal of Computer Assisted Tomography 2(1) (1978), 100–108.

Nalcioglu

, Lou

R.Y.

Post-reconstruction method for beam hardening in computerised tomography, Physics in Medicine and Biology 24(2) (1979), 330–340.

Tang

, Mou

, Xu

, et al., Data consistency condition–based beam-hardening correction, Optical Engineering 50(7) (2011), 076501.

Tang

, Mou

, Yang

et al., Application of projection simulation based on physical imaging model to the evaluation of beam hardening corrections in X-ray transmission tomography, Journal of X-ray Science and Technology 16(2) (2008), 95–117.

10.

Tang

, Zhang

, Wang

, et al., An optimal bone correction capable of simultaneously suppressing various types of image artifacts, CT Theory and Applications 27(3) (2018), 301–314.

11.

Chen

, Zhang

, Xia

, et al., Non-convex primal-dual algorithm for image reconstruction in spectral CT, Computerized Medical Imaging and Graphics 87 (2021), 101821.

12.

Lauritsch

, Bruder

Head phantom, (1998) [2018]. http://www.imp.uni-erlangen.de/phantoms/head/head.html.

13.

Fewell

T.R.

, Shuping

R.E.

, Healy

K.E.

Handbook of Computed Tomography X-ray Spectra. Washington, DC: HHS Publication (FDA), 1981, 8162:27–29.

14.

Rowlands

J.A.

, Taylor

K.W.

Absorption and noise in cesium iodide X-ray image intensifiers, Medical Physics 10(6) (1983), 786–795.

15.

Wang

, Lin

, Cheng

, et al., A general cone-beam reconstruction algorithm, IEEE Transaction on Medical Imaging 12(3) (1993), 486–496.

16.

Kalender

W.A.

, Hebel

, Ebersberger

Reduction of CT artifacts caused by metallic implants, Radiology 164(2) (1987), 576–577.

17.

Hsieh

Adaptive streak artifact reduction in computed tomography resulting from excessive X-ray photon noise, Medical Physics 25(11) (1998), 2139–2147.

18.

Liao

, Lin

W.A.

, Zhou

S.K.

, et al., ADN: artifact disentanglement network for unsupervised metal artifact reduction, IEEE Transactions on Medical Imaging 39(3) (2019), 634–643.

19.

Zhang

, Yu

Convolutional neural network based metal artifact reduction in x-ray computed tomography, IEEE Transactions on Medical Imaging 37(6) (2018), 1370–1381.

20.

Basu

, Breseler

Uniqueness of tomography with unknown view angles, IEEE Transactions on Image Processing 9(6) (2000), 1094–1106.

21.

Natterer

The mathematics of computerized tomography, Philadelphia: SIAM, 2001.

22.

, Mou

, Cai

Calibration of cone beam rotational X-ray image sequence, Chinese Journal of Electronics 13(3) (2004), 424–428.