Sparse-view X-ray CT based on a box-constrained nonlinear weighted anisotropic TV regularization

1.   Introduction

In clinical applications, X-ray computed tomography (CT) aims to visualize the internal structure of the human body by reconstructing the tissues' attenuation coefficients $\mu$ to X-rays. Depending on diverse X-ray sources, there are parallel beam, fan beam, and cone beam CTs ^[1,2,3,4]. In this paper, for ease of explanation, we focus on image reconstructions in parallel beam CT, even though the proposed method can be used in other fan beam and cone beam CTs. In a parallel beam CT, parallel X-ray beams in different directions are transmitted through the patient, who lies between the X-ray sources and the detectors (see Figure 1). The corresponding attenuated X-ray intensities are measured through the detectors. The inverse problem of a parallel beam CT is to reconstruct the attenuation coefficient from the received attenuated X-ray intensities. Given the perfect measured data, the reconstruction methods include filtered back-projection, the algebraic reconstruction technique (ART) ^[1], and so on. Due to the advantages of speed, accuracy, and excellence with bones and lungs for nondestructive testing, X-ray CT is widely used in medical imaging to aid doctors in diagnosing diseases ^[1].

Nevertheless, exposure of the patients to the environment of radiation increases the risk of many diseases such as leukemia, cancer, etc ^[5]. Low dose CT is an effective way to reduce such a risk. There are generally three methods of low dose CT. The first is to reduce the tube voltage/currents, the second is the limited angle CT reconstruction, and the third is the sparse-view CT reconstruction. For the first method, to obtain a meaningful CT image, we need an efficient denoising algorithm ^[6]. For the second method, visible singularities, invisible singularities, and artifacts exist^[1]. For the third method, since it violates the Nyquist/Shannon sampling criterion, strong streaking artifacts will occur^[7]. In this paper, we focus on removing the streaking artifacts in the sparse-view CT reconstruction. To this end, we need to develop efficient ways to deal with the ill-posedness caused by the projection downsamplings ^[1].

To deal with the ill-posedness, the data-driven and model-driven methods exist. The data-driven method includes the usage of a convolutional neural network (CNN) ^[8], U-Net ^[9], and the GoogLeNet ^[10]. However, as pointed out in ^[11], there should be more evidence of such methods being used in clinical applications. For the sparse-view CT, it is also very difficult to provide enough labeled data, since we cannot produce the data with full projections without enough doses. Hence, the model-driven method is still quite necessary.

For the model-driven method, note that the algebraic method ^[12] has the flexibility of incorporating the a-priori information of $\mu$. It is widely used in a sparse-view CT. To be precise, the reconstruction of $\mu$ is recast into a minimization problem of an energy functional constructed by a (weighted) least square fitting and a regularization.

Noting the piecewise constant structure of medical images, its gradient can be considered sparse. The study performed in ^[13] proposed the total variation (TV) regularization method in sparse-view CT. However, it is well known that TV will introduce new blocky/staircasing artifacts ^[14]. The study performed in ^[15] proposed the anisotropic TV regularization method, that could produce distortions along axes. To handle such problems, many variations of TV regularization have been proposed in the last two decades. The study performed in ^[16] proposed an edge-preserving TV regularizer which used the $e^{-|\nabla \mu |^2/\sigma^2}$ as the edge detector, where $\sigma$ was a prescribed parameter that represented the amount of smoothing. Later, ^[17] proposed a similar discretized version to ^[16]. The similar methods can be found in ^[18,19]. Nevertheless, the ability of removing the streaking artifact and edge-preserving can be further improved due to the fact that the amount of regularization near the edges and away from that does not differ much since $e^{-|\nabla \mu |^2/\sigma^2}\in [0, 1]$. The study performed in ^[20] proposed the total generalized variation (TGV) regularizer, which used the second order derivative of the unknown $\mu$. While this method can avoid the blocky artifacts, it assumes the piecewise linear structure of the image. It is well known that the piecewise constant is commonly used for a medical image. The study performed in ^[21] proposed a directional TV regularizer in which a directional derivative was considered, rather than just using ${\nabla} \mu$. However, different weights can employed in different directions to improve the performance of this method. The study performed in ^[22] proposed the $L^p (0 < p < 1)$ regularization method. However, the images heavily relied on the parameter $p$.

$L_{1}/L_{2}$ is a recently proposed regularization technique ^[23]. This method is based on updating the regularization parameter in each iteration. However, the updating parameter is not region-dependent, that is, in each iteration, the minimization problem is isotropic. The study performed in ^[24] proposed a nonlinear weighted anisotropic TV (NWATV) regularization method and used it in electrical impedance tomography, which is a low resolution imaging modality. In this paper, a box-constrained NWATV method was used in sparse-view CT, which produced a significantly improved reconstruction compared with directly using NWATV and box-constrained $L_{1}/L_{2}$ methods. Precisely, across the internal edges where ${\nabla}\mu \to \infty$, we set the regularization to be small to preserve the edge, while near the smooth region, we set a normal regularization to make the ill-posed problem better posed. We found a significant convergence behavior of the iteration process with the box constraint (set the reconstruction value lies in a proper interval). We validated the proposed algorithm using the Shepp-Logan phantom, the walnut X-ray data provided by Finnish Inverse Problems Society (http://fips.fi/dataset.php), and the clinical lung image provided by The Cancer Imaging Archive (TCIA: https://www.cancerimagingarchive.net/collection/lungct-diagnosis/).

The rest of the paper is organized as follows. In Section 2, we provide a brief introduction of the parallel beam CT. In Section 3, we introduce the proposed box-constrained nonlinear weighted anisotropic TV regularization and provide an iterative reconstruction algorithm. In Section 4, we validate the performance of the proposed regularization method using the Shepp-Logan phantom, the actual walnut CT experiment data, and the clinical lung image. In Section 5, we discuss the rules of the choice of regularization parameters. In Section 6, we conclude the paper and provide some future research topics.

2.   Preliminaries of parallel beam CT

In parallel beam CT, we restrict our explanations to the two-dimensional space because the parallel beam always lies in a plane, and each time, the X-ray can only pass through a slice of the object. Let $\Omega\subset\mathbb{R}^2$ represent a bounded region of the imaging object. Denote $\mu$ as the attenuation coefficient of $\Omega$, which is generally a piecewise constant function in medical imaging. In parallel beam CT, an incident X-ray beam along the direct lines $L_{\theta, s} : = \{\mathbf{x}\in \mathbb{R}^2\; :\; \Theta \cdot \mathbf{x} = s\}$ passes through the object, which lies between the X-ray sources and the detectors (see Figure 1). Here, $s\in \mathbb{R}$ denotes the signed distance of $L_{\theta, s}$ to the original point $O(0, 0)$ and $\Theta = (\cos \theta, \sin \theta)$, with $\theta\in [0, \pi)$ denoting the angle of a directly line $l$ and $x$-axis, where $l$ is perpendicular to $L_{\theta, s}$. We assume that the incident X-ray intensity is $I^0$. For a fixed $\theta$ and $s$, the attenuated X-ray intensity $I(\theta, s)$ can be measured through the detector. The relation between the measured $I(\theta, s)$ and the unknown $\mu$ is described by the Lambert-Beer law ^[1,2]:

where $\mathcal{R}_\theta[f](s)$ is the Radon transform ^[25] of $f$, which is defined as follows:

with $d\ell_{\mathbf{x}}$ denoting the length element.

In medical imaging, we assume that a parallel beam contains $J$ X-rays, and hence $J$-detectors are employed to detect the corresponding attenuated X-rays. We assume that the $J$ X-rays are equidistantly distributed. To be precise, we assume that the signed distances of the X-rays to the original point lie in $[\underline{s}, \bar s]$, and the signed distance of the $j$-th $(j\in \{1, 2, \cdots, J\})$ X-ray to the original point is $s_j = \underline s + (j-1) (\overline s - \underline s) /(J-1)$. Moreover, we utilize $K$ angles $\theta_k = (k-1)\pi/K$ for $1 \leq k \leq K$ and $k\in \mathbb{Z}_+$. Then, through the detectors $M: = JK$, projection datum $I(\theta_k, s_j)$ ($1\leq j\leq J$ and $1\leq k\leq K$) can be measured.

For ease of explanation, we denote $y_m = \ln\frac{I^0}{I(\theta_k, s_j)}$ for $m = J(k-1) + j$. Then, we have the following:

for $1\leq m\leq M$. Suppose the field of view ($\mbox{FOV}$) is the square, and $\mbox{FOV} = [a, b]\times [a, b]$ which satisfies $\Omega \subset \mbox{FOV}$. We discretize $\mbox{FOV}$ to be $N\times N$ pixels $P_{qt}$ ($1\leq q, t \leq N$). Then, from Eq (2.1), the reconstruction of $\mu$ can be recast to solve the following linear system:

Here, $A = (a_{mp})$ is an $M \times n$ matrix for $n = N^2$, $\mathbf{y} = (y_m)$ is an $M\times 1$ vector, and $\mathbf{u} = (u_p)$ is an $n\times 1$ vector. To be precise, $a_{mp}$ is the length of the projection line which lies in the pixel $P_{qt}$, i.e., $a_{mp} = |L_{\theta_k, s_j}\cap P_{qt}|$ for $p = N(t-1)+q$, $u_p = \mu(q, t)$ and $y_m$ is defined in Eq (2.1).

3.   Nonlinear weighted anisotropic TV regularization with box constraint

Note that for the sparse-view CT, we generally have $M\ll n$; hence, to solve Eq (2.2), we reformulate it to the following least squares problem:

where $\|\cdot\|_{\ell_{2}}$ denotes the standard Euclidean norm in $\mathbb{R}^{M}$. Since $A^TA$ is ill-conditioned, where $\cdot^T$ represents the transpose of $\cdot$, we approximate Eq (3.1) by the following well-conditioned problem:

The first term on the right-hand side of Eq (3.2) is the data fidelity term, the second term ${\mbox{Reg}}(\mathbf{u})$ is the regularization term, and $\lambda$ is the regularization parameter that balances the fidelity and the regularization terms. Precisely, instead of seeking the solution in the space $\ell_2$, where there may be infinitely many solutions, we seek the solution in its subspace, which is characterized by ${\mbox{Reg}}(\mathbf{u})$.

Since the edge of the internal structure is a key feature in medical imaging, the choice of the term ${\mbox{Reg}}(\mathbf{u})$ should obey the following rule:

● Near the local edges of $\mu$ or $\mathbf{u}$ where $|{\nabla} \mu|\approx \infty$, do as little regularization as possible to preserve the edges.

● Perform a normal regularization when $\mu$ is smooth (i.e., $|{\nabla}\mu|$ is small to make Eq (3.2) well-posed).

● The range of the reconstructed $\mathbf{u}$ coincides with the range of its true value.

Combining the above three considerations, we define ${\mbox{Reg}}(\mathbf{u})$ as follows:

where $\gamma$ is the indicator of using box constraint, that is, $\gamma = 1$ if the box constraint is used, while $\gamma = 0$ if no box constraint is employed. Here, $\mathbf{p} = \left(\omega(\mathbf{\mathcal{D}}_{x}\mathbf{u}); \omega(\mathbf{\mathcal{D}}_{y}\mathbf{u})\right)\in\mathbb{R}^{2n}, \mathbf{\mathcal{D}} \mathbf{u} = (\mathbf{\mathcal{D}}_{x}\mathbf{u}; \mathbf{\mathcal{D}}_{y}\mathbf{u})\in\mathbb{R}^{2n}$, where $\mathbf{\mathcal{D}}_{x}$, $\mathbf{\mathcal{D}}_{y}\in\mathbb{R}^{n\times n}$ are the first-order difference operators along the $x$ and $y$ directions, respectively, $\omega(\cdot) = \frac{1}{|\cdot|^{2}+\beta}$ with $\beta > 0$ a small positive number to avoid zero being the denominator; $\Pi_{[c_1, c_2]}(\mathbf{u})$ is an indicator function, which equals to $0$ if for all $i\in \{1, 2, \cdots, n\}$, where $\mathbf{u}[i]\in [c_1, c_2]$, and equals to $+\infty$ otherwise. Note that $\Pi_{[c_1, c_2]}$ is capable of enforcing $\mathbf{u}$ into the range of the actual attenuation coefficient.

The augmented Lagrangian functional of Eq (3.2), together with Eq (3.3), can be expressed as follows:

where $\mathbf{d}, \mathbf{v}$ are the auxiliary variables, $\mathbf{b}, \mathbf{e}$ are the Lagrangian multipliers, and $\rho, \alpha$ are the scalar penalty parameters.

To minimize Eq (3.4), we use the alternating direction method of multipliers (ADMM)^[26]. To be precise, for an initial guess ($\mathbf{d}^{(0)}, \mathbf{p}^{(0)}, \mathbf{v}^{(0)}, \mathbf{b}^{(0)}, \mathbf{e}^{(0)}$), $\mathbf{u}$ is iteratively updated via the following scheme:

For Eq (3.5a) and Eq (3.5b), the minimizers have the following closed form:

and

Here, $\mathcal{I}$ is an $n\times n$ identity matrix; $\mathbf{b}^{(k)}[i]$, $\mathbf{p}^{(k)}[i]$, $\mathbf{u}^{(k+1)}[i]$, and $\mathbf{d}^{(k+1)}[i]$ are the $i$-th element of $\mathbf{b}^{(k)}$, $\mathbf{p}^{(k)}$, $\mathbf{u}^{(k+1)}$, and $\mathbf{d}^{(k+1)}$, respectively, and $h_{g}(\cdot)$ represents the soft threshold formula, which is defined as follows^[27]:

where $\mathrm{sgn}$ is the sign function.

We end this section by summarizing the aforementioned process as the reconstruction algorithm in the form of the pseudocode shown in Algorithm 1.

4.   Experiments

In this section, to validate the advantages of the proposed regularization, we perform experiments using the Shepp-Logan numerical model, the walnut actual CT model provided by Finnish Inverse Problems Society (http://fips.fi/dataset.php), and the clinical CT image provided by The Cancer Imaging Archive (https://www.cancerimagingarchive.net/collection/lungct-diagnosis/). Additionally, the walnut CT data is provided in ZENODO (https://zenodo.org/record/1254206). The experimental models are shown in Figure 2.

4.1.   Experiment setup

To show the advantages of the proposed regularization, we compare the reconstructions of the most recently proposed gradient-based $L_1/L_2$^[23] with box constraint and the nonlinear weighted anisotropic TV regularization^[24] with box constraint. To compare the performance of the reconstructions, we compute the relative errors (including the $L^2$ relative errors ${\mbox{RE}}(k)$ and the $H^1$ relative errors $\widetilde{{\mbox{RE}}}(k)$) and mean square errors ${\mbox{MSE}}(k)$ for the $k$-th step as follows:

and

Here, $\mathbf{u}^{(k)}$ represents the result of the reconstruction at the $k$-th step. In the Shepp-Logan phantom, $\mathbf{u}_0$ represents the ground truth image; in the walnut experiment, it represents the CT image reconstructed using full projections and the filtered back-projection (FBP) method, since we do not know the ground truth image in the actual experiment; in the lung image, it represents the actual image from TCIA. Figure 2(a)–(c) illustrate $\mathbf{u}_0$ of the Shepp-Logan phantom, the walnut, and the lung experiments.

Besides, we also compare the peak signal-to-noise ratio ${\mbox{PSNR}}(k)$ and the structural similarity index ${\mbox{SSIM}}(k)$ ^[28] for the $k$-th step, which is defined as follows:

and

where $\cdot\odot\cdot$ represents the componentwise multiplication, $\mu_k$ and $\sigma_k$ represent the local mean and the local standard deviation of $\mathbf{u}^k$, respectively, and $\mu_0$ and $\sigma_0$ represent the local mean and the local standard deviation of $\mathbf{u}_0$, respectively; moreover, $\sigma_{k0}$ denotes the cross-covariance of $\mathbf{u}^{(k)}$ and $\mathbf{u}_{0}$, and $C_{1} = 10^{-4}, C_{2} = 9\times 10^{-4}$ are set to be the default values in the Matlab build-in function "ssim".

In the numerical experiment, we set the size of the reconstructed images to be $256\times 256$ pixels and set the number of detectors to be $J = 362$. In the reconstruction of the actual walnut experiment, the size of the reconstructed images is set to be $164 \times 164$, and the number of detectors is also 164. In the lung image, we use the first image of patient R_172 in the dataset. The original size is 512 $\times$ 512, though we evenly sample to get the ground truth image of 128 $\times$ 128. The number of detectors is set to be 181. Using the Radon transform, we obtain the projection data corresponding to the lung image. To solve Eq (3.6), we use the generalized minimal residual algorithm (GMRES) ^[29] to accelerate the computation.

The reconstructions are carried out using Matlab 2018a (The MathWorks, Inc., Natick, MA, USA) on a workstation with 1.60 GHz Inter (R) Core (TM) i5-8250U CPU, 8.00 GB memory, Windows 10 operating system. Additionally, we use the MATLAB package AIR Tools Ⅱ to simulate the parallel beam for the CT scanning ^[30].

4.2.   Numerical experiment results

Under box constraints, we compare the results of $L_{1}/L_{2}$ method and nonlinear weighted TV regularization. In all experiments, we set the maximum number of external and internal iterations in the box-constrained $L_{1}/L_{2}$ to be 300 and 5, respectively. For fair comparisons, the ranges of other parameters are set according to ^[23] to minimize the $L^{2}$ relative error (RE) and to obtain the best performance. The number of iterations in NWATV(without box constraint) and the box-constrained NWATV regularization is set to be 300. We set the candidate set of parameters to be $\lambda\in\left\{0.002, 0.004, 0.006, 0.008, 0.01\right\}$, $\rho\in\{20, 40, 60,200,400,600\}$, and $\alpha\in\left\{5, 20, 40, 60\right\}$. The selection of parameters has been carefully optimized to achieve a balance between minimizing the RE and optimizing the visual effects. Henceforth, in the figures and tables, we use $L_1/L_2$-box and NWATV-box as the abbreviations of the box-constrained $L_{1}/L_{2}$ and the box-constrained NWATV, respectively.

First, we consider the effect of box constraint for the NWATV regularization. We consider a parallel beam CT reconstruction with 31 angles uniformly taken from 0º to 150º, and the noise level is 0.5%. The box constraint is [0, 1]. Therefore, the sample size is set to be 362 $\times$ 31. We do our best to choose the parameters such that RE attains the minimum value. Precisely, we choose $\rho$ = 20, $\lambda$ = 0.002, $\alpha$ = 5 in the box-constrained NWATV, and $\rho$ = 20, $\lambda$ = 0.004 in NWATV. The reconstruction images are shown in Figure 3, and the corresponding relative errors (RE and $\widetilde{{\mbox{RE}}}$), mean square errors (MSE), peak signal-to-noise rations (PSNR), SSIM values, and the CPU times are shown in Table 1.

The results show that with a short reconstruction time, the box-constrained NWATV regularization can perform better than that without the box constraint. In Figure 4, we illustrate the evolutions of ${\mbox{RE}}(k)$ and ${\mbox{SSIM}}(k)$ with $k$, the iteration step. The figure clearly shows that the box constraint improves the convergence behavior of the NWATV method.

Next, we show that the proposed regularization can reconstruct a satisfied image using different sampling angles and is robust against the Gaussian random noise. We evenly take 90, 60, and 30 angles for comparisons from 0º to 179º. For each case, we add different Gaussian random noises with the levels 0.5%, 1%, 1.5%, and 2%. The box constraint is set to be [0, 1].

We list the specific parameters selected in Table 2. The corresponding reconstruction results are shown in Figures 5–7. We show the values of the corresponding numerical results in Table 3.

As we can see from the reconstructions, the box-constrained NWATV can achieve numerical results similar to that of the box-constrained $L_{1}/L_{2}$ method at each sampling size. However, the CPU time is reduced at least 8 times. Taking visual effects into account, in the case of small angles (362 $\times$ 30), as the noise level increases, the box-constrained NWATV has more advantages than the box-constrained $L_{1}/L_{2}$ in noise removal and detail recovery.

4.3.   Walnut CT experiment and human lung experiment

First, in this section, we validate the performance of the proposed regularization method using the walnut data from the experimental datasets^[31], which includes its complete projection matrix $\bar{A}\in\mathbb{R}^{19680\times 26896}$ and the projection data $\bar{\mathbf{y}}\in\mathbb{R}^{164\times 120}$ (3º per projection). We evenly subsample 50 angles from 0º to 149º by employing the first 50 projections. Then, the corresponding projection matrix $A$ is also subsampled such that $A\in\mathbb{R}^{8200\times 26896}$. We use a high-resolution filtered back-projection reconstruction^[31], and then evenly sample the obtained image to obtain the reference image $\mathbf{u}_{0}$, which is shown in Figure 2(b). From the reference image, we can estimate the box constraint is [0, 0.6]. Considering the magnitudes of $A^{T}A$ and $\mathbf{\mathcal{D}}^{T}\mathbf{\mathcal{D}}$, the parameters we take here are listed as follows: $\rho = 0.4$, $\lambda = 0.0001$, and $\alpha = 5$. For the parameters in the box-constrained $L_{1}/L_{2}$, we balance noise removal and detail recovery, choosing a set of parameters with the best visual effects. The results are shown in Figure 8, and the values of the corresponding numerical results are shown in Table 4. The arrows in Figure 8 depict that the box-constrained NWATV method perfoms better than the box-constrained $L_1/L_2$ in edges and details preserving.

Finally, we discuss the application of the two methods to clinical data^[32]. We uniformly sample the image of 512 $\times$ 512 to obtain the image of 128 $\times$ 128, which is shown in Figure 2(c), and we represent it as the column vector $\mathbf{u}_0\in\mathbb{R}^{128^2}$. We consider parallel beam scanning and set the sampling sizes as 181 $\times$ 60 and 181 $\times$ 30 to obtain the projection matrix A. The projection is $\mathbf{y} = A\mathbf{u}_0$. Since there is already noise in the lung image, we do not need to add additional noise any more. The box constraint is taken as [0, 0.06]. In the experiment of box-constrained NWATV, we set the parameters as $\lambda = 0.001$, $\rho = 3\times 10^{-9}$, and $\alpha = 0.01$ for both the 30 and 60 degrees projections. As in the box-constrained $L_1/L_2$, we balance noise removal and detail recovery, choosing a set of parameters to make the reconstructed image have the best visual effect. The results are shown in Figure 9, and the corresponding numerical results are shown in Table 5.

For further comparison, we also illustrate profiles of the reconstructed walnut and human lung images along the dash lines shown in Figure 10. From the profiles, it's clear that the proposed model behaves better in edges and details preserving.

In conclusion, through the aforementioned experiments, the box-constrained NWATV method can produce similar and visually better results than the recently developed box-constrained $L_1/L_2$ method, while the CPU time is significantly decreased.

5.   Discussion

Sparse-view CT and limited angle CT are two important ways to reduce the risk of radiation exposure in medical CT scanning. The recent development of artificial intelligence promotes the medical applications of low dose CT ^[33]. However, recent research has reported the instabilities of such methods ^[11,34]. Moreover, to gather the high quality training data, the conventional regularization-based reconstruction methods are still quite necessary. This is because to get the labeled data, without using an elegant regularization, data obtained from the full projections should be employed, which exposes the patient under high risk of radiation.

Another possible model-based approach is the sinogram inpainting method ^[35]. However, it will cause other artifacts. Hence, proper regularization in the image reconstruction process is a mild way to produce high performance CT images for diagnoses and to gather the training data for AI methodology. ^[36,37] provide other ways of combining the model-based and data-driven methods for CT image reconstructions.

5.1.   Parameter selection rules

In regularization-based reconstructions, the choice of the regularization parameter is generally very difficult and criticized. However, from the point view of dimensional reduction, we can produce a high performance CT image with a much higher dimension given a small amount of parameters. On the one hand, for the choices of the parameters, mathematically there are the discrepancy principle methods and the statistical methods ^[38] to deal with proper choices of the parameters. On the other hand, there are some rules for the parameters, which are listed as follows.

● As in ^[24], the performance of the reconstruction (i.e., the relative error depends only the $\frac{\lambda}{\rho}$ rather than the values of $\lambda$ and $\rho$).

● The values of $\rho$ and $\alpha$ should be determined by order of $A^{T}A$ and the number of pixels in such a way that $A^{T}A+\rho\mathbf{\mathcal{D}}^{T}\mathbf{\mathcal{D}}+\alpha\mathcal{I}$ should not change its order.

● The optimal $\frac{\lambda}{\rho}$ value ranges from $10^{-5}$ to $10^{-4}$, since we can minimize the relative error by using such a value.

● The optimal range of $\lambda$ is approximately $10^{-3}$ to $10^{-2}$.

In Figure 11, we depict the evolution of RE with $\frac{\lambda}{\rho}$ and $\alpha$ for different $\rho$. From the figures, we can see that for each $\rho$, RE is invariant with the changes of $\alpha$ when $\frac{\lambda}{\rho}$ is fixed.

5.2.   Comparisons on the computation cost with $L_1/L_2$ method

In this section, we explain the computation load of the proposed box-constrained NWATV method and compare with that in box-constrained $L_1/L_2$ method ^[23]. To guarantee the convergence of the box-constrained $L_1/L_2$ method, both inner and outer loops have to be used, while in the proposed NWATV method, a single loop could produce a good performance of the convergence. Note that in each loop, the most expensive computation is the calculation in (3.6). Similar to ^[23], in this paper, the maximum number of inner loops is set to be $5$, which means that for the same outer loops, the computational load is at least 5 times bigger than that in the proposed method.

6.   Conclusions and future works

In this paper, we proposed a box-constrained nonlinear weighted anisotropic TV regularization method and used it in sparse-view CT. Using the Shepp-Logan phantom and the actual walnut models, we validated that the proposed regularization could reconstruct a more accurate CT image than the most recently developed $L_{1}/L_{2}$ regularization method. To be precise, the reconstruction time was reduced more than 8 times while maintaining similar relative errors and a structural similarity index. The proposed method showed advantages, especially when the sampling angles were less than 60 and the noise level was more than 1%. Additionally, numerical simulations displayed a good convergence performance of the proposed iterative scheme. Moreover, since the pixel values of digital images were mostly limited to a certain range, it was reasonable to add box constraints in image processing ^[1,23]. Note that errors/noises in the iteration based numerical scheme may accumulate with the iterations, and the box constraint plays the role in suppressing the accumulation in some extent. Hence, the box constraint could enforce the iteration converges to the critical point of the functional $\mathcal{L}$, as shown in Figure 4. We note that in the box-constrained $L_1/L_2$ method, the authors also note the similar role of the box constraint ^[23].

Future works should definitely include the mathematical theory of the convergence of the proposed iterative scheme. Moreover, in this paper, we list some rules for the selection of parameters $\lambda$, $\rho$, and $\alpha$ through manual tuning. In the future, we could develop an automatic way to select the optimal parameters by minimizing the discrepancy function in an admissible set using the combinatorial optimization method. If the noise level $\delta$ is known, the discrepancy function could be defined as $F(\lambda, \rho, \alpha) = |\|A\mathbf{u}^{\lambda, \rho, \alpha} - \mathbf{y}\|-\tau \delta|$ for a given $\tau > 1$ to avoid underregularization. On the other hand, if the noise level is not known, the Hanke Raus function $F(\lambda, \rho, \alpha) = \|A\mathbf{u}^{\lambda, \rho, \alpha} - \mathbf{y}\|/\lambda$ could be used ^[39]. To minimize the above functions, we first need to select a good initial guess $\lambda_0, \rho_0, \alpha_0$, and the optimal parameters could be obtained by minimizing the discrepancy functional using an alternate direction iteration scheme. Furthermore, the proposed method can be further used in the area of metal artifact reduction (MAR) in CT reconstruction ^[40], beam-hardening artifact reduction ^[41], limited angle artifact reduction ^[1], and so on.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by Shandong Provincial Outstanding Youth Fund (Grand No. ZR2018JL002), National Natural Science Foundation of China (Grand No. 12271312).

Conflict of interest

The authors declare there is no conflict of interest.