Box-Constrained Nonlinear Weighted Anisotropic TV Regularization for Beam Hardening Artifacts Reduction in CT

Xue Shi

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/jamp.2025.132019                

Abstract

In Computed Tomography (CT), the beam hardening artifacts are caused by polychromatic X-ray beams applied in real medical imaging. In this article, we applied the recently proposed box-constrained nonlinear weighted anisotropic total variation regularization (box-constrained NWATV) method in the process of the reconstruction. We do numerical experiments to validate the advantages of the proposed method in reducing the beam hardening artifacts compared with the existing ways.

Keywords

Beam Hardening Artifacts, Box-Constrained NWATV, Computed Tomography

1. Introduction

Because of the change in voltage and current, the energy of incident X-rays is changed, which makes polychromatic beams in CT scanning. Since low-energy photons are more easily absorbed than high-energy photons, this causes the beam to harden as X-rays pass through the object [1]. The incident intensity is denoted by $I^0\left(E\right)$ and $I_{\theta }\left(E,s\right)$ is the measured attenuated intensity along $L_{\theta ,s}$ . So, the average energy of the X-rays reaching the detectors $\frac{{\displaystyle \int E{I}_{\theta }\left(E,s\right)\text{d}E}}{{\displaystyle \int I_{\theta }\left(E,s\right)\text{d}E}}$ is higher than that of the incident X-rays $\frac{{\displaystyle \int E{I}^0\left(E\right)\text{d}E}}{{\displaystyle \int I^0\left(E\right)\text{d}E}}$ . Then, from reference [2], we acquire the inequality

$\frac{{\displaystyle \int E{I}^0\left(E\right)\text{d}E}}{{\displaystyle \int I^0\left(E\right)\text{d}E}}\le \frac{{\displaystyle \int E{I}_{\theta }\left(E,s\right)\text{d}E}}{{\displaystyle \int I_{\theta }\left(E,s\right)\text{d}E}}$

This effect is known as beam hardening artifacts. The effects will cause the cupping artifacts in the reconstructed attenuation images.

Beam hardening artifact reduction is important in medical and industrial CT applications to improve the visual quality of images [3]. Since 1975, the correction of beam hardening artifacts has been a hot topic in CT. The correction methods are divided into two categories: hardware correction methods and software correction methods. The hardware correction methods add some correction tools to the CT system to suppress the beam hardening artifacts. Software correction methods are based on the mathematical view of beam hardening artifacts. Common software correction methods mainly include the polynomial fitting method [4], Monte Carlo correction method [5], dual-energy method [6], iterative correction method [7], and single-energy correction method [3].

Iterative reconstruction method models the problem mathematically and transforms the model into an optimization problem with fidelity term and regularization term. In reference [8], the authors proposed the Nonlinear Weighted Anisotropic TV (NWATV) regularization in the electrical impedance tomography to solve the EIT inverse problem. In reference [9], the authors proposed the box-constrained nonlinear weighted anisotropic TV (box-constrained NWATV) regularization to solve the sparse-view X-ray CT inverse problem.

In this article, we proposed an iterative method to correct the beam hardening artifacts. We build a discretized model of beam hardening artifacts and use the box-constrained NWATV regularization to correct beam hardening artifacts and compare the reconstructed results of box-constrained NWATV regularization, TV regularization, and ISP method by some numerical experiments to validate the advantages of the box-constrained NWATV regularization.

2. Notation and Concept

In CT imaging, the imaging object, such as the human body, is placed between the X-ray sources and the detectors. The X-rays are injected into the body and the attenuated X-rays are measured on the detector. X-ray imaging visualizes the internal structure of the object by reconstructing the attenuation coefficients via the relationship between the injection and measurements described by the law of Lambert-Beer.

Suppose the object is located in a two-dimensional bounded region Ω, and for simplicity, we assume the parallel beams are used in this paper. To be precise, $\theta$ distribute evenly from 0 to $\frac{179\pi }{180}$ . Suppose 180 angles are used. We discretize $\text{Ω}$ to be $N\times N$ pixels. We assume that there are $K$ rays at each angle and $K$ detectors to receive the signals. We assume the sign distance between each X-ray and the origin of the region Ω is $\left[{\mathcal{S}}_1,{\mathcal{S}}_2\right]$ , and the distance from the origin of the kth ray in angle $\theta$ is $s_k={\mathcal{S}}_1+\left(k-1\right)\cdot \frac{\left({\mathcal{S}}_2-{\mathcal{S}}_1\right)}{K-1}$ . Then, the number of all X-rays is $M=180\cdot K$ . So, the matrix of parallel beam scanning $A$ is $M\times N^2$ , and the measured data $y$ is an $M\times 1$ vector.

The kth monochromatic X-ray beam in angle $\theta$ passes through the object along $L_{\theta ,s_k}$ and the law of Lambert-beer is described as [10]

$I_{\theta ,s_k}=I^0\cdot {\text{e}}^{-{\displaystyle {\int }_{L_{\theta ,s_k}}\mu \text{d}l}}.$

When X-ray beam passes through the object along $L_{\theta ,s_k}$ from the low-energy $E_{\text{min}}$ to high-energy $E_{\text{max}}$ , the incident intensity is $I^0={\displaystyle {\int }_{E_{\text{min}}}^{E_{\text{max}}}{I}^0\left(E\right)\text{d}E}$ , and $\mu \left(E\right)$ is the attenuation coefficient of the object [11]. Then, the relationship between the outgoing intensity and the attenuation coefficient is described by the law of Lambert-Beer [11]

$I_{\theta ,s_k}={\displaystyle {\int }_{E_{\text{min}}}^{E_{\text{max}}}{I}^0\left(E\right)\cdot \text{exp}\left\{-{\displaystyle {\int }_{L_{\theta ,s_k}}\mu \left(E\right)\text{d}l}\right\}\text{d}E}.$

The measured data of the mth X-ray is

$y_m=-\text{ln}\left(\frac{I_{\theta ,s_k}}{I^0}\right)=-\text{ln}\left({\displaystyle {\int }_{E_{\text{min}}}^{E_{\text{max}}}\frac{I^0\left(E\right)}{{\displaystyle {\int }_{E_{\text{min}}}^{E_{\text{max}}}{I}^0\left(E\right)\text{d}E}}\text{exp}\left\{-{\displaystyle {\int }_{L_{\theta ,s_k}}\mu \left(E\right)\text{d}l}\right\}\text{d}E}\right),$

where $m=\theta \cdot K+k$ .

3. Polychromatic Model and Algorithm

We build the Lambert-Beer law for polychromatic X-rays as a discretized mathematical model. We discretize the energy $E=\left[E_1,E_2,\cdots ,E_J\right]$ of incident X-rays and get the equation ${\displaystyle {\int }_{E_{\text{min}}}^{E_{\text{max}}}{I}^0\left(E\right)\text{d}E}={\displaystyle \sum }_{i=1}^J{I}^0\left(E_i\right)$ . So, we obtain the measured data of the mth X-ray

$y_m=-\text{ln}\left({\displaystyle \sum }_{i=1}^J\frac{I^0\left(E_i\right)}{{{\displaystyle \sum }}_{i=1}^J{I}^0\left(E_i\right)}\text{exp}\left\{-{\displaystyle {\int }_{L_{\theta ,s_k}}\mu \left(E_i\right)\text{d}l}\right\}\right)$

and

${\text{e}}^{-y_m}={\displaystyle \sum }_{i=1}^J\frac{I^0\left(E_i\right)}{{{\displaystyle \sum }}_{i=1}^J{I}^0\left(E_i\right)}\text{exp}\left\{-{\displaystyle {\int }_{L_{\theta ,s_k}}\mu \left(E_i\right)\text{d}l}\right\}.$

The right-hand side of the equation is equivalent to a weighted average, resulting in a monochrome image with the energy of $E_0$ with $E_{\text{min}}<E_0<E_{\text{max}}$ . Therefore, the above equation is equivalent to $y_m={\displaystyle {\int }_{L_{\theta ,s_k}}\mu \left(E_0\right)\text{d}l}$ , that is, the attenuation coefficient of the object to be scanned at this time is the attenuation coefficient when the energy $E_0$ is attenuated. From this, we get the following reconstruction model $y=Au\left(E_0\right)$ , where $u\left(E_0\right)$ is the image, we need to reconstruct, i.e.

$u^{\text{*}}\left(E_0\right)=\text{arg}\underset{u\left(E_0\right)}{\text{min}}{\Vert y-Au\left(E_0\right)\Vert }_{l_2}^2.$

In reference [9], we get the nonlinear weighted anisotropic TV regularization and the equation

$u^{\text{*}}\left(E_0\right)=\text{arg}\underset{u\left(E_0\right)}{\text{min}}\frac{1}{2}{\Vert y-Au\left(E_0\right)\Vert }_{l_2}^2+\lambda {\Vert p\cdot \mathcal{D}u\left(E_0\right)\Vert }_{l_1}+\frac{\gamma }{\lambda }{\displaystyle \prod }_{\left[c_1,c_2\right]}\left(u\left(E_0\right)\right),$

where $\lambda$ is the regularization parameter. ${\displaystyle {\prod }_{\left[c_1,c_2\right]}\left(u\left(E_0\right)\right)}$ is an indicator function. It will be 0 if $c_1\le u\left(E_0\right)\left[i\right]\le c_2$ , and will be $+\infty$ otherwise. If $\gamma =1$ means box constraint is used otherwise $\gamma =0$ . And $p=\left(\omega \left({\mathcal{D}}_{x}u\left(E_0\right)\right);\omega \left({\mathcal{D}}_{y}u\left(E_0\right)\right)\right)\in {ℝ}^{2N^2}$ , $\mathcal{D}u\left(E_0\right)=\left({\mathcal{D}}_{x}u\left(E_0\right);{\mathcal{D}}_{y}u\left(E_0\right)\right)\in {ℝ}^{2N^2}$ , $\omega (\cdot )=\frac{1}{{\left|\cdot \right|}^2+\beta }$ , where the $\beta$ is a positive number to avoid ${\left|\cdot \right|}^2$ to be 0. The augmented Lagrangian functional is

$\begin{array}{l}L\left(u\left(E_0\right),d,p,v;e,b\right)\\ =\frac{1}{2}{\Vert Au\left(E_0\right)-y\Vert }_{l_2}^2+\lambda {\Vert p\cdot d\Vert }_{l_1}+\langle b,\mathcal{D}u\left(E_0\right)-d\rangle \\ +\frac{\rho }{2}{\Vert \mathcal{D}u\left(E_0\right)-d\Vert }_{l_2}^2+{\displaystyle \prod }_{\left[c_1,c_2\right]}\left(v\right)+\langle e,u\left(E_0\right)-v\rangle +\frac{\alpha }{2}{\Vert u\left(E_0\right)-v\Vert }_{l_2}^2\end{array}$ (1)

where the $e,b$ are the Lagrangian multipliers and the $\alpha ,\rho$ are the penalty parameters [9]. The Alternative Direction Multiplier Method (ADMM) [12] is used to minimize (1). To be precise, we update $u\left(E_0\right)$ , $d$ , $p$ , $b$ , $v$ and $e$ by (2)-(8). They are updated as

$u{\left(E_0\right)}^{\left(n+1\right)}={\left(A^{\text{T}}A+\alpha I+\rho {\mathcal{D}}^{\text{T}}\mathcal{D}\right)}^{-1}\left(A^{\text{T}}y+{\mathcal{D}}^{\text{T}}{b}^{\left(n\right)}-\rho {\mathcal{D}}^{\text{T}}{d}^{\left(n\right)}+e^{\left(n\right)}-\alpha v^{\left(n\right)}\right),$ (2)

$d^{\left(n+1\right)}\left[i\right]=h_{\frac{\lambda \left|p^{\left(n\right)}\left[i\right]\right|}{\rho }}\left(\mathcal{D}u{\left(E_0\right)}^{\left(n+1\right)}\left[i\right]+\frac{1}{\rho }{b}^n\left[i\right]\right)$ (3)

where $b^{\left(n\right)}\left[i\right]$ , $p^{\left(n\right)}\left[i\right]$ , $u^{\left(n+1\right)}\left[i\right]$ , $d^{\left(n+1\right)}\left[i\right]$ represent the $i\text{th}$ element of $b^{\left(n\right)}$ , $p^{\left(n\right)}$ , $u^{\left(n+1\right)}$ , $d^{\left(n+1\right)}$ . The $h_g$ is defined as

$h_g(\cdot )=\{\begin{array}{ll}\cdot -g\mathrm{sgn}(\cdot )\hfill & \left|\cdot \right|>g;\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$ (4)

And the $p$ , $b$ , $v$ and $e$ are defined as

$p^{\left(n+1\right)}=\left(\omega \left({\mathcal{D}}_{x}u{\left(E_0\right)}^{\left(n+1\right)}\right);\omega \left({\mathcal{D}}_{y}u{\left(E_0\right)}^{\left(n+1\right)}\right)\right),$ (5)

$b^{\left(n+1\right)}=b^n+\rho \left(\mathcal{D}u{\left(E_0\right)}^{n+1}-d^{\left(n+1\right)}\right),$ (6)

$v^{\left(n+1\right)}=\text{min}\left\{\text{max}\left\{u{\left(E_0\right)}^{\left(n+1\right)}+\frac{1}{\alpha }{e}^{\left(n\right)},c_1\right\},c_2\right\},$ (7)

$e^{\left(n+1\right)}=e^{\left(n\right)}+\alpha \left(u{\left(E_0\right)}^{\left(n+1\right)}-v^{\left(n+1\right)}\right)$ (8)

with the initial data $d^{\left(0\right)}=0$ , $p^{\left(0\right)}=\left(\frac{1}{\beta }\right)1$ , $v^{\left(0\right)}=0$ , $e^{\left(0\right)}=0$ , $b^{\left(0\right)}=0$ . We set

the maximum iteration to 1000, unless the iterative is broken by $\Vert u^{\left(n+1\right)}-u^{\left(n\right)}\Vert <\epsilon$ . And we get the matrix $A$ of pallel beam scanning in CT by the MATLAB package AIR Tools II [13].

4. Experiments

To display the advantages of the box-constrained NWATV regularization, we will compare the reconstructed results via the Peak Signal-to-Noise Ratio (PSNR), and Structural Similarity Index (SSIM) of box-constrained NWATV regularization, TV regularization in reference [14] and ISP method in reference [1]. We use the MATLAB built-in function “SSIM” to get the results of SSIM. $\text{PSNR}\left(n\right)$ is defined as

$\text{PSNR}\left(n\right)=10{\log }_{10}\frac{\text{max}\left(u_n\cdot \odot \cdot u_n\right)}{\text{MSE}\left(n\right)},$

where $\cdot \odot \cdot$ means component-wise multiplication [8]. The $u^{\text{*}}$ means the ground truth image, $u_n$ represents the reconstruction in the nth iteration, and

$\text{MSE}\left(n\right)=\frac{1}{N^2}\left({\Vert u_n-u^{\text{*}}\Vert }_{l_2}^2\right).$

In numerical experiments, we use the 128 × 128 standard Shepp-Logan phantom and the 100 × 100 three-disk model made by Matlab in Figure 1 as the ground truth image for box-constrained NWATV regularization and TV regularization. Because of the speed of the calculation, we use the downsampled 32 × 32 standard Shepp-Logan phantom and 50 × 50 three-disk model as the ground truth image for ISP method. In three-disk model, we get the reconstructed images and the $\text{PSNR}\left(n\right)$ , $\text{SSIM}\left(n\right)$ by box-constrained NWATV regularization with parameters $\alpha =100$ , $\beta ={10}^{-10}$ , $\lambda =0.2$ , $\rho =60$ , TV regularization and the ISP method in Figure 2 and Table 1. From Figure 2, we can see that the box-constrained NWATV regularization is more complete for the preservation of boundary information.

Table 1. From left to right are the PSNR and SSIM for box-constrained NWATV regularization, TV regularization, and ISP method for Figure 2.

NWATV-box		TV		ISP
PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
29.2330	0.6895	32.1228	0.5639	14.4854	0.1480

Figure 1. The left side is the 128 × 128 Shepp-Logan ground truth image and the right is the 100 × 100 three-disk model.

Figure 2. The first row of images from left to right are the image with beam hardening artifacts and the reconstructed images for box-constrained NWATV regularization, the TV regularization, and the ISP method. Images in the second row are profiles of reconstructed images for box-constrained NWATV regularization and TV regularization in the 25th, 50th, and 80th lines, and ISP method in the 13th, 25th, and 40th lines.

In Shepp-Logan phantom numerical experiments, we obtain the $\text{PSNR}\left(n\right)$ , $\text{SSIM}\left(n\right)$ by box-constrained NWATV regularization with $\alpha =100$ , $\beta =1\times {10}^{-10}$ , $\lambda =0.002$ , $\rho =800$ , TV regularization and the ISP method in Figure 3 and Table 2. From the results of SSIM, it can be seen that the results of the reconstruction of the box-constrained NWATV regularization are closer to the ground truth images than the other two methods.

Table 2. From the left to right are the PSNR and SSIM for box-constrained NWATV regularization, TV regularization and ISP method for Figure 3.

NWATV-box		NWATV		ISP
PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
11.0017	0.7024	11.7279	0.6578	23.7687	0.0959

Figure 3. The meaning of each figure in the first row is the same as that in Figure 2. The images in the second row from the left to right are profiles of reconstructed images for box-constrained NWATV regularization and TV regularization in the 35th lines and ISP method in the 9th lines.

5. Conclusion

From the above results, the NWATV-box regularization has obvious advantages in eliminating the beam hardening artifacts. According to the reconstructed images, we know that TV regularization will blur the image and to a certain extent, will weaken the details, which is not the case in the approach we propose. Moreover, because the ISP method requires the gradient descent method, so the box-constrained NWATV regularization has a higher speed of calculation than ISP method. We can see that the data of SSIM in box-constrained NWATV regularization are better than ISP method and TV regularization.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Van Gompel, G., Van Slambrouck, K., Defrise, M., Batenburg, K.J., de Mey, J., Sijbers, J., et al. (2011) Iterative Correction of Beam Hardening Artifacts in CT. Medical Physics, 38, S36-S49. https://doi.org/10.1118/1.3577758
[2]	Seo, J.K. and Woo, E.J. (2012) Nonlinear Inverse Problems in Imaging. Wiley https://doi.org/10.1002/9781118478141
[3]	Nalcioglu, O. and Lou, R.Y. (1979) Post-Reconstruction Method for Beam Hardening in Computerised Tomography. Physics in Medicine & Biology, 24, 330-341. https://doi.org/10.1088/0031-9155/24/2/009
[4]	Herman, G.T. (1979) Correction for Beam Hardening in Computed Tomography. Physics in Medicine & Biology, 24, 81-106. https://doi.org/10.1088/0031-9155/24/1/008
[5]	Gang, Z. (2006) Beam Hardening Correction Based on Monte Carlo Simulation. High Energy Physics and Nuclear Physics, 30, 178-182. https://api.semanticscholar.org/CorpusID:138325575
[6]	Marshall, W.H., Alvarez, R.E. and Macovski, A. (1981) Initial Results with Prereconstruction Dual-Energy Computed Tomography (Predect). Radiology, 140, 421-430. https://doi.org/10.1148/radiology.140.2.7255718
[7]	Elbakri, I.A. and Fessler, J.A. (2002) Statistical Image Reconstruction for Polyenergetic X-Ray Computed Tomography. IEEE Transactions on Medical Imaging, 21, 89-99. https://doi.org/10.1109/42.993128
[8]	Song, Y., Wang, Y. and Liu, D. (2022) A Nonlinear Weighted Anisotropic Total Variation Regularization for Electrical Impedance Tomography. IEEE Transactions on Instrumentation and Measurement, 71, 1-13. https://doi.org/10.1109/tim.2022.3220288
[9]	Li, H. and Song, Y. (2024) Sparse-View X-Ray CT Based on a Box-Constrained Nonlinear Weighted Anisotropic TV Regularization. Mathematical Biosciences and Engineering, 21, 5047-5067. https://doi.org/10.3934/mbe.2024223
[10]	Brabant, L., Pauwels, E., Dierick, M., Van Loo, D., Boone, M.A. and Van Hoorebeke, L. (2012) A Novel Beam Hardening Correction Method Requiring No Prior Knowledge, Incorporated in an Iterative Reconstruction Algorithm. NDT & E International, 51, 68-73. https://doi.org/10.1016/j.ndteint.2012.07.002
[11]	Hansen, P.C., Jørgensen, J.S. and Lionheart, W.R.B. (2021) Computed Tomography: Algorithms, Insight, and Just Enough Theory. Society for Industrial and Applied Mathematics. https://epubs.siam.org/isbn/978-1-61197-666-3
[12]	Boyd, S. (2010) Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in Machine Learning, 3, 1-122.
[13]	Hansen, P.C. and Jørgensen, J.S. (2017) AIR Tools II: Algebraic Iterative Reconstruction Methods, Improved Implementation. Numerical Algorithms, 79, 107-137. https://doi.org/10.1007/s11075-017-0430-x
[14]	Rudin, L.I., Osher, S. and Fatemi, E. (1992) Nonlinear Total Variation Based Noise Removal Algorithms. Physica D: Nonlinear Phenomena, 60, 259-268. https://doi.org/10.1016/0167-2789(92)90242-f