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

: Sparse-view computed tomography (CT) is an important way to reduce the negative e ff ect of radiation exposure in medical imaging by skipping some X-ray projections. However, due to violating the Nyquist / Shannon sampling criterion, there are severe streaking artifacts in the reconstructed CT images that could mislead diagnosis. Noting the ill-posedness nature of the corresponding inverse problem in a sparse-view CT, minimizing an energy functional composed by an image fidelity term together with properly chosen regularization terms is widely used to reconstruct a medical meaningful attenuation image. In this paper, we propose a regularization, called the box-constrained nonlinear weighted anisotropic total variation (box-constrained NWATV), and minimize the regularization term accompanying the least square fitting using an alternative direction method of multipliers (ADMM) type method. The proposed method is validated through the Shepp-Logan phantom model, alongisde the actual walnut X-ray projections provided by Finnish Inverse Problems Society and the human lung images. The experimental results show that the reconstruction speed of the proposed method is significantly accelerated compared to the existing L 1 / L 2 regularization method. Precisely, the central processing unit (CPU) time is reduced more than 8 times.


Introduction
X-ray computed tomography (CT) aims to visualize the internal structure of the human body by reconstructing tissues' attenuation coefficients µ to X-rays in clinical applications.Depending on diverse X-ray sources, there are parallel beam, fan beam and cone beam CT [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 CT.In 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 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 fast, accurate and excellence with bones and lungs for nondestructive testing, X-ray CT is widely used in medical imaging to aid doctors diagnose diseases [1].
Nevertheless, the exposure of 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 ways 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 sparseview CT reconstruction.For the first one, to obtain a meaningful CT image, we need an efficient denoising algorithm [6].For the second one, there exist visible singularities, invisible singularities and artifacts [1].For the third one, 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 datadriven method includes the usage of 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 can not 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 µ.It is widely used in sparse-view CT.To be precise, the reconstruction of µ is recast into a minimization problem of an energy functional constructed by a (weighted) least square fitting and a regularization.
Noting that the piecewise constant structure of medical images, its gradient can be considered sparse.[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].[15] proposed the anisotropic TV regularization method while it could produce distortions along axes.To handle such problems, many variations of TV regularization have been proposed in the last two decades.[16] proposed an edge-preserving TV regularizer which used the e −|∇µ| 2 /σ 2 as the edge detector, where σ is a prescribed parameter representing 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 improved more due to the fact that the amount of regularization near the edges and away from that does not differ much since e −|∇µ| 2 /σ 2 ∈ [0, 1].[20] proposed the total generalized variation (TGV) regularizer which takes use of the second order derivative of the unknown µ.While this method can avoid the blocky artifacts, it assumes the piecewise linear structure of the image.It is well known that commonly for a medical image it should be piecewise constant.[21] proposed a directional TV regularizer in which a directional derivative is considered rather than just using ∇µ.However, different weights can employed in different directions to improve the performance of this method.[22] proposed L p (0 < p < 1) regularization method.However, the images are 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.[24] proposed a nonlinear weighted anisotropic TV (NWATV) regularization method and used it in electrical impedance tomography, a low resolution imaging modality.In this paper, a box-constrained NWATV method is used in sparse-view CT which produces a significantly improved reconstruction compared with directly using NWATV and box-constrained L 1 /L 2 methods.Precisely, across the internal edges where ∇µ → ∞, 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 validate 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/lungctdiagnosis/).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.

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 Ω ⊂ R 2 represent a bounded region of the imaging object.Denote µ as the attenuation coefficient of Ω, which is generally a piecewise constant function in medical imaging.In parallel beam CT, an incident X-ray beam along the direct lines L θ,s := {x ∈ R 2 : Θ • x = s} passes through the object which lies between the X-ray sources and detectors (see Figure 1).Here s ∈ R denotes the signed distance of L θ,s to the original point O(0, 0) and Θ = (cos θ, sin θ) with θ ∈ [0, π) denoting the angle of a directly line l and x-axis, where l is perpendicular to L θ,s .We assume that the incident X-ray intensity is I 0 .For fixed θ and s, the attenuated X-ray intensity I(θ, s) can be measured through the detector.The relation between the measured I(θ, s) and the unknown µ is described by the Lambert-Beer law [1,2]  where R θ [ f ](s) is the Radon transform [25] of f defined as with dℓ 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 [s, s] and the signed distance of the j-th For ease of explanation, we denote y m = ln I 0 I(θ k ,s j ) for m = J(k − 1) + j.Then we have Then from the Eq (2.1), the reconstruction of µ can be recast to solve the following linear system

.2)
Here A = (a mp ) is an M × n matrix for n = N 2 , y = (y m ) is an M × 1 vector and u = (u p ) is an n × 1 vector.To be precise, a mp is the length of the projection line which lies in the pixel P qt , i.e.
+ q, u p = µ(q, t) and y m is defined in Eq (2.1).

Nonlinear weighted anisotropic TV regularization with box constraint
Note that for the sparse-view CT, generally we have M ≪ n, hence to solve the Equation (2.2), we reformulate it to the following least squares problem where ∥ • ∥ ℓ 2 denotes the standard Euclidean norm in R M .Since A T A is ill-conditioned, where • T represents the transpose of •, 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 Reg(u) is the regularization term and λ is the regularization parameter which balances the fidelity and the regularization terms.Precisely, instead of seeking the solution in the space ℓ 2 where there may be infinitely many solutions, we seek the solution in its subspace characterized by Reg(u).
Since the edge of internal structure is a key feature in medical imaging, the choice of the term Reg(u) should obey the following rule: • Near the local edges of µ or u where |∇µ| ≈ ∞, do as little regularization as possible to preserve the edges.
• The range of the reconstructed u coincides with the range of its true value.
Combining the above three considerations, we define Reg(u) as follows where γ is the indicator of using box constraint, that is γ = 1 if the box constraint is used while , where D x , D y ∈ R n×n are respectively the first-order difference operators along the x and y with β > 0 a small positive number to avoid zero being the denominator; and equals to +∞ otherwise.Note that Π [c 1 ,c 2 ] is capable of enforcing 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 where d, v are the auxiliary variables, b, e are the Lagrangian multipliers and ρ, α are the scalar penalty parameters.
We end this section by summarizing the above process as the reconstruction algorithm in the form of pseudocode shown in Algorithm 1.

Experiments
In this section, to validate the advantages of the proposed regularization we do experiments The walnut CT data is also provided in ZENODO (https://zenodo.org/record/1254206).The experimental models are shown in Figure 2 .

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 Algorithm 1 The box-constrained NWATV method Require: Projection matrix A, observed data y, and a bound [c 1 , c 2 ] for the original image.
end if 12: end for 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 RE(k) and the H 1 relative errors RE(k)) and mean square errors MSE(k) for the k-th step as follows Here, u (k) represents the result of the reconstruction at k-th step.In Shepp-Logan phantom, Besides we also compare the peak signal-to-noise ratio PSNR(k) and the structural similarity index SSIM(k) [28] for the k-th step defined as follows In the numerical experiment, we set the size of the reconstructed images to be 256 × 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 × 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×512, we sample evenly to get the ground truth image of 128×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 Equation (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.We also make use of the MATLAB package AIR Tools II to simulate parallel beam for the CT scanning [30].

Numerical experiment results
Under box constraints, we compare the results of L 1 /L 2 method and nonlinear weighted TV regularization.In all the 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]  We first consider the effect of box constraint for the NWATV regularization.We consider 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×31.We do our best to choose the parameters such that RE attains the minimum value.Precisely, we choose ρ=20, λ=0.002, α=5 in the box-constrained NWATV, and ρ=20, λ=0.004 in NWATV.
The reconstruction images are shown in Figure 3, and the corresponding relative errors (RE and 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 RE(k) and 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 for 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 the Table 2.The corresponding reconstruction results are shown in Figure 5, 6, 7. We show the values of the corresponding numerical results in Table 3.   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.

Mathematical Biosciences and Engineering
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 above 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.[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 expose 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, a proper regularization in the image reconstruction process is a mild way to produce high performance CT images for diagnosis and 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.

Parameter selection rules
In regularization-based reconstructions, the choice of 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 by given a small amount of parameters.For the choices of the parameters, on the one hand, mathematically there are the methods as discrepancy principle and 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 λ ρ rather than the values of λ and ρ.
• The values of ρ and α should be determined by order of A T A and the number of pixels in such a way that A T A + ρD T D + αI should not change its order.
• The optimal λ ρ value ranges from 10 −5 to 10 −4 since using such value we can minimize the relative error.
In Figure 11 we depict the evolution of RE with λ ρ and α for different ρ.From the figures, we can see that for each ρ, RE is invariant with the changes of α when λ ρ is fixed.

Comparisons on the computation cost with
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).
As that in [23], the maximum number of inner loop is set to be 5 in this paper which means that for the same outer loops, the computational load is at least 5 times bigger than that in the proposed method.The range of λ is {0.002, 0.004, 0.006, 0.008, 0.01, 0.02, 0.04, 0.06, 0.08, 0.1}.models, we validate that the proposed regularization can reconstruct a more accurate CT image than the most recently developed L 1 /L 2 regularization method.To be precise, the reconstruction time is reduced more than 8 times while maintaining similar relative errors and structural similarity index.The proposed method shows advantages, especially when the sampling angles are less than 60 and the noise level is more than 1%.Numerical simulations also show a good convergence performance of the proposed iterative scheme.Moreover, since the pixel values of digital images are mostly limited to a certain range, it is 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, 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 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 through manual tuning, we list some rules for the selection of parameters λ, ρ, α.In the future, we could develop an automatic way of the selection of optimal parameters by minimizing the discrepancy function in an admissible set using the combinatorial optimization method.If the noise level δ is known, the discrepancy function could be defined as F(λ, ρ, α) = |∥Au λ,ρ,α − y∥ − τδ| for a given τ > 1 to avoid underregularization.

Mathematical Biosciences and Engineering
Volume xx, Issue x, xxx-xxx y∥/λ could be used [39].To minimize the above functions, we need firstly select a good initial guess λ 0 , ρ 0 , α 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.

MathematicalFigure 1 .
Figure 1.Schematic diagram of a parallel X-ray beam CT system.
0 represents the ground truth image; in walnut experiment, it represents the CT image reconstructed using full projections and the method of filtered back-projection (FBP) since we do not know the ground truth image in the actual experiment; in lung image, it represents the actual image from TCIA.Figure2 (a)-(c) illustrate u 0 of Shepp-Logan phantom, walnut and lung experiment.

Figure 2 .
Figure 2. Experimental models we used.(a) is the Shepp-Logan phantom model in numerical experiment, (b) shows the walnut model in actual CT experiment which is obtained in Finnish Inverse Problems Society (FINNISH), and (c) is a clinical lung image provided by the Cancer Image Archive (TCIA).

Figure 3 .
Figure 3.Effect of box constraint in NWATV reconstruction.The grayscale window is [0,1].Left: Reconstruction result under box constraint.Right: Reconstruction result without box constraint.

Figure 4 .
Figure 4.The effect of box constraint on the convergence of the iteration scheme of NWATV.The left is the RE and the right is the SSIM.

MathematicalFigure 5 .Figure 6 .
Figure 5. Reconstruction results using the box-constrained NWATV and the boxconstrained L 1 /L 2 methods with the sampling size 362×90.The top row is the results of the box-constrained L 1 /L 2 , while the bottom row is the results of the box-constrained NWATV.From left to right are respectively the reconstructed results with noise levels of 0.5%, 1%, 1.5% and 2%.

Figure 9 .
Figure 9. Reconstruction results of the human lung image.The first row depicts the reconstructed images using the box-constrained L 1 /L 2 (left) and box-constrained NWATV (right) methods, respectively for the sampling size of 181×60.The second row illustrates the similar results for sampling size of 181×30.

Figure 10 .
Figure 10.The profiles of the reconstructed walnut and human lung images along the dash lines on left figures.In the right figures, the red line illustrate the values of u 0 along the dash line, the blue line show the values of reconstructions using box-constrained L 1 /L 2 methods while the green line depict the reconstructions using box-constrained NWATV methods along the dash lines.The sampling size of the human lung imaging is set to be 181×30.
to minimize the L 2 relative error(RE) and obtain the best performance.The number of iterations in NWATV(without box constraint) and Mathematical Biosciences and Engineering Volume xx, Issue x, xxx-xxx the box-

Table 1 .
Numerical results of the Logan-Shepp model using the NWATV reconstruction with and without box constraint.
angles(362×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.

Table 3 .
Performance of the box-constrained L 1 /L 2 and the box-constrained NWATV regularization for Shepp-Logan phantom for different sampling sizes and noise levels.

Table 4 .
Comparisons of the performances of two box-constrained methods using the walnut data.

Table 5 .
Comparisons of the performances of two box-constrained methods using the human lung image.