 Research
 Open Access
 Published:
Sparse CT reconstruction based on multidirection anisotropic total variation (MDATV)
BioMedical Engineering OnLine volume 13, Article number: 92 (2014)
Abstract
Background
The sparse CT (Computed Tomography), inspired by compressed sensing, means to introduce a prior information of image sparsity into CT reconstruction to reduce the input projections so as to reduce the potential threat of incremental Xray dose to patients’ health. Recently, many remarkable works were concentrated on the sparse CT reconstruction from sparse (limitedangle or fewview style) projections. In this paper we would like to incorporate more prior information into the sparse CT reconstruction for improvement of performance. It is known decades ago that the given projection directions can provide information about the directions of edges in the restored CT image. ATV (Anisotropic Total Variation), a TV (Total Variation) norm based regularization, could use the prior information of image sparsity and edge direction simultaneously. But ATV can only represent the edge information in few directions and lose much prior information of image edges in other directions.
Methods
To sufficiently use the prior information of edge directions, a novel MDATV (MultiDirection Anisotropic Total Variation) is proposed. In this paper we introduce the 2DIGS (Two Dimensional Image Gradient Space), and combined the coordinate rotation transform with 2DIGS to represent edge information in multiple directions. Then by incorporating this multidirection representation into ATV norm we get the MDATV regularization. To solve the optimization problem based on the MDATV regularization, a novel ART (algebraic reconstruction technique) + MDATV scheme is outlined. And NESTA (NESTerov’s Algorithm) is proposed to replace GD (Gradient Descent) for minimizing the TVbased regularization.
Results
The numerical and real data experiments demonstrate that MDATV based iterative reconstruction improved the quality of restored image. NESTA is more suitable than GD for minimization of TVbased regularization.
Conclusions
MDATV regularization can sufficiently use the prior information of image sparsity and edge information simultaneously. By incorporating more prior information, MDATV based approach could reconstruct the image more exactly.
Background
CT (Computed Tomography) is one of the most important medical image technologies. Due to the potential cancer risk associated with the radiation dose in CT, recent technical development focuses on reducing radiation dose in CT. Several CT machine manufacturers have developed CT machines with iterative reconstruction algorithms to reduce radiation dose, such as the IRIS (Iterative Reconstruction in Image Space) of Siemens, ASiR (Adaptive Statistical Iterative Reconstruction) of GE (General Electric Co.), iDose (a trademark) of Philips and AIDR (Adaptive Iterative Dose Reduction) 3D of Toshiba, which are reported to reduce the radiation dose by 3050% while maintaining the reconstructed image quality comparable to the conventional FBP (Filtered Back Projection) method [1–5].
On the other hand, biomedical engineering researchers are trying to improve the performance of the iterative reconstruction algorithms by introducing more prior information of the reconstructed images to reduce the projection input and the radiation dose even more. One important prior information is image sparsity, whose usefulness has been testified by the recent popular compressed sensing technique in signal processing [6–8]. Before the compressed sensing, the sparsity of medical images has been used for limitedangle tomography [9] and fewview CT bloodvessel reconstruction [10]. Then following the compressed sensing theory, two kinds of sparse models were proposed for the sparse CT reconstruction. By ‘sparse CT’ we mean two aspects. First, a single CT image is approximately sparse like other natural images, and the information changes between successive CT images in dynamic CT are also approximately sparse. Second, the projection data is sparsely collected, such as limitedangle tomography [9] and fewview CT [10], rather than the full view dense collection in conventional CT scan. The first model uses the sparsity coming from the gradient transformation of the image which is assumed to be approximately piecewiseconstant, such as ARTTVMIN (Algebraic Reconstruction Technique Total Variation MINimization) [11] and ASDPOCS (Adaptive Steepest Descent Projection Onto Convex Sets) [12]. The second model uses the sparsity coming from the subtraction of the reconstructed image from its prior image, such as PICCS (Prior Image Constrained Compressed Sensing) [13], PIRPLE (Prior Image Registered Penalized Likelihood Estimation) [14, 15] and FCCS (Feature Constrained Compressed Sensing) [16].
In this paper we consider the methods based on the first model for the sparse CT reconstruction problems. In this model, the reconstruction problem can be written as a constrained optimization
The equality constraint is a guarantee of the data fidelity, where ℳ is the projection matrix modeling the forward projection, $\overrightarrow{f}$ is the image vector to be restored, and $\overrightarrow{g}$ is the projection vector. The objective function is a norm function of $\overrightarrow{f}$, which is a regularization for introducing some prior information like image sparsity. By using norm function of $\overrightarrow{f}$, we want to obtain a sparse solution of this constrained optimization. Therefore, sparser $\overrightarrow{f}$ is, smaller ${\overrightarrow{f}}_{\mathit{\text{reg}}}$ should be. A commonly used norm function is TV (Total Variation) norm [11]. The optimization problem can also be written as an unconstrained form
where $\lambda $ is a regularization parameter adjusting the relative penalties on the sparseness and the data fidelity.
To promote the performance of sparse CT reconstruction, two kinds of improvements are often considered. The first kind of improvement focuses on the sparsity regularization. Sidky [17] proposed the sparser TpV (Total pVariation) norm to replace the TV norm. Yu [18] used the Haar wavelet transformation to make use of sparsity. Jia [19] replaced the TV norm with tight frame and used GPU (Graphics Processing Unit) to accelerate the reconstruction speed. Apart from replacement of TV norm, Tian [20], Liu [21], and Chang [22] developed some adaptive adjustment factors embedded in the TV norm to enhance the edge characteristic. Recently, Liu [23] introduced the TVStokes regularizer into sparse CT, which can smooth the isophote in the target image so as to preserve the detail and smoothness of the reconstructed image. The second kind of improvement focuses on the algorithms for resolving (2), which include CG (Conjugate Gradient) [24], forward backward splitting plus GPU [25], splitbregman [26], GPBB (Gradient Projection Barzilai Borwein) plus GPU [27], ChamboleePock algorithm [28], and ADMM (Alternating Direction Method of Multipliers) [29]. Besides, Niu and Zhu [30] proposed to use the logbarrier term to approximate the data fidelity term, which removes the forward and backward projection in the traditional iterative algorithms like ART (Algebraic Reconstruction Technique).
In this paper, we consider the first kind of improvement. The MDATV (MultiDirection Anisotropic Total Variation) is proposed to replace TV norm. The aim of MDATV is to introduce into sparse CT reconstruction more prior information of edges, as the ATV (Anisotropic Total Variation) does [31, 32]. But ATV can only represent the edge information in few directions, which loses much prior information of edge directions. This paper combines the coordinate rotation transform with 2DIGS (Two Dimensional Image Gradient Space), and represents the edge information in multiple directions easily. Based on this multidirection representation, MDATV can incorporate more edge information into the sparse CT reconstruction, so as to improve the performance.
The rest of the paper is organized as follows. In section Methods, firstly, the background of ART and ATV is reviewed. Then the MDATV norm and its corresponding minimization methods are described. Finally, the ART + MDATV scheme is outlined. Section Simulations and experiments uses the proposed MDATV based approach to solve sparse CT reconstruction problems in numerical simulations and with experimental data, so as to demonstrate its remarkable efficiency. Finally, we conclude with section Discussions and conclusions, discussing the pros and cons of MDATV.
Methods
Review of ART + ATV
ART
The ART + ATV is proposed based on ART + TV method [11, 12]. ART updates the estimated image iteratively. Firstly, the estimated image is forward projected into the sinogram space, then the difference between the estimated sinogram and the given projections is back projected into the image space to update the estimated image. These steps are operated iteratively until some termination criteria are met. This method is also known as POCS (Projection Onto Convex Sets) in linear algebra. The ART update formula is as following
where g _{ m } is the m th element of the projection $\overrightarrow{g}$. ${\overrightarrow{M}}_{m}$ is the m th row vector of the projection matrix ℳ. $\overrightarrow{f}\left[n,m\right]$ is the estimated image after fusing the m th element of the projection data during the n th iteration. The superscript T represents transpose.
TV and ATV
The TV norm in the ART + TV method is used to introduce and strengthen the prior information of image sparsity. This prior information significantly improves the quality of reconstructed image with low computation load. The TV norm is defined as
where f _{ i,j } is the pixel on the i th row and the j th column in the target image. D _{ i } f _{ i,j } =f _{ i,j } f _{ i–1,j }, and D _{ j } f _{ i,j } =f _{ i,j } f _{ i,j–1 } are the vertical and the horizontal gradients respectively, D _{ i } and D _{ j } are the corresponding finite differential operators. Because the image edges are orthogonal to their corresponding gradients, D _{ i } f _{ i,j } and D _{ j } f _{ i,j } can also represent edges in the horizontal and the vertical directions respectively. The TV norm is isotropic, since it assigns the same energy to both the vertical and the horizontal gradients. Different from TV norm, ATV norm is anisotropic and defined as
where A and B are weights controlling the energy ratios of the vertical and the horizontal gradients, respectively. When A = B, it is isotropic, else it is anisotropic. (4) can also be written as
where η=A/B is the weight factor. Therefore, we only need one parameter to control the energy ratios of the vertical and the horizontal gradients without affecting the results of optimization. In this paper we set η=1000. Based on experiments, this value of η can make the edges in the enhanced direction strong enough, while the potential artifacts in the suppressed direction weak enough. When η is smaller, the suppression to the potential artifacts may be not sufficient. A bigger η is also valid, but the reconstructions hardly change any more.
Motivation of ATV
The motivation of ATV is to introduce into sparse CT reconstruction a prior information of edge directions. This prior information was discovered by Quinto [33] decades ago. In sparse CT reconstruction, the edge information tangent to the projection rays would be measured and restored easily, while those not tangent to the projection rays should be harder to ‘see’. This phenomenon can be theoretically explained by the central slice theory [34].
Without loss of generality, we take the parallel projection for example. This example is shown in Figure 1. Suppose the projection rays in direction a produce a row of projections. Each point of the projection row corresponds to one projection ray, and each point’s value is the linear integration of the attenuation coefficients passed through by the corresponding projection ray. According to the central slice theory, the 1D Fourier transformation of this projection row is the central slice perpendicular to a in the Fourier frequency space. In the 2D Fourier space, the central slice is a line passing through the origin. Thereby the projection data collected by projection rays in direction a corresponds to a central slice perpendicular to a in the 2D Fourier space.
According to the Fourier optics [35], if a family of edges in a 2D image are in the direction a, then the Fourier transformation of these edges is a line passing through the origin and perpendicular to a in the 2D Fourier frequency space, as shown in Figure 2. This means that the line passing through the origin and perpendicular to a in the 2D Fourier space corresponds to the edge information in direction a of the 2D image.
Compare Figures 1 and 2, we find that the projection rays in direction a record the edge information in the same direction. Therefore, the edge information tangent to the projection rays would be restored easily, since these edge information is well recorded by the projection rays. While the edge information perpendicular to the projection rays is not recorded, thereby these edges are hard to ‘see’.
The discussion above indicates that the geometrical information of projection may affect the sparse CT reconstruction. For example in the limitedangle projection, if most projections are tangent or approximately tangent to the horizontal direction, then the isotropic TV regularization may result in some artifacts in the vertical direction, which will blur the clear edges in the horizontal direction. Since the projections did not record any edge information in the vertical direction, the energy vacancy in the vertical direction would be filled by artifacts. However, ATV could handle this problem by assigning more energy in the horizontal direction and less energy in the vertical direction, thus to enhance the edge information in the horizontal direction and suppress the artifacts in the vertical direction.
MDATV
However, ATV can only represent edges in few directions. This may lose much prior information of edges in many other directions. To remedy this defect, MDATV is proposed to use the entire prior information of edges.
2DIGS
Before description of MDATV, we firstly introduce the 2DIGS. In a 2D discrete image, f _{ i,j } is defined as the pixel on the i th row and the j th column. The 2DIGS $\mathbf{\Psi}$ contains two parts. One is the vertical edges denoted by the 2D matrix E _{ v }, and the other is the horizontal edges denoted by the 2D matrix E _{ h }. They are defined as
where D _{ v } and D _{ h } are the same as D _{ i } and D _{ j } in (3) representing the vertical and the horizontal finite differential operators respectively. The E _{ h } and E _{ v } of SheppLogan phantom are shown in Figure 3(b) and (e). The original image is generated by the function ‘phantom(128)’ in MATLAB^{®}, as shown in Figure 3(a).
Edges in multiple directions
E _{ h } and E _{ v } represent the horizontal and the vertical edges of the 2D image, which are used in the TV and ATV norm. To represent edges in many more other directions, the coordinate rotation transform is introduced to be applied to the 2DIGS.
First, we recall the elementary coordinate rotation transform in 2D space. Figure 4 shows the transform process. The coordinate system xOy rotates counterclockwise by θ degrees to the coordinate system x′Oy′, where 0^{°} ≤ θ ≤ 180^{°}. The target point T in xOy is denoted by (x _{0}, y _{0}). Its corresponding rotated target point T′ is denoted by $\left({x}_{0}^{\prime},{y}_{0}^{\prime}\right)$ in x′Oy′. Obviously, we have ${x}_{0}={x}_{0}^{\prime}$ and ${y}_{0}={y}_{0}^{\prime}$. When we need to process the rotated target point T′ in xOy, we can project the coordinates of T′ in x′Oy′ onto the coordinate system xOy. Let (p, q) denote the projected coordinates of T′ in xOy. According to the elementary geometry, the transform formula is
where p is the linear combination of the projections of ${x}_{0}^{\prime}$ and ${y}_{0}^{\prime}$ onto the x axis. q is the linear combination of the projections of ${x}_{0}^{\prime}$ and ${y}_{0}^{\prime}$ onto the y axis.
When the coordinate rotation transform is applied to the 2DIGS $\mathbf{\Psi}$, the horizontal edges E _{ h } corresponds to the horizontal line segment x _{0} in (6), and the vertical edges E _{ v } corresponds to the vertical line segment y _{0} in (6). After counterclockwise rotation by $\theta $ degrees, we get the rotated 2DIGS Ψ′, which contains the rotated horizontal and vertical edges ${\mathit{E}}_{\mathit{h}}^{\mathit{\prime}}$ and ${\mathit{E}}_{\mathit{v}}^{\mathit{\prime}}$. But we need to process the rotated horizontal and vertical edges in the 2DIGS $\mathbf{\Psi}$. Therefore, we need to project ${\mathit{E}}_{\mathit{h}}^{\mathit{\prime}}$ and ${\mathit{E}}_{\mathit{v}}^{\mathit{\prime}}$ onto the 2DIGS $\mathbf{\Psi}$. Let E _{ hr } and E _{ vr } denote the projected ${\mathit{E}}_{\mathit{h}}^{\mathit{\prime}}$ and ${\mathit{E}}_{\mathit{v}}^{\mathit{\prime}}$. Then E _{ hr } and E _{ vr } correspond to p and q in (6). Thus we have
By adjusting the angle $\theta $, the edges in any direction could be denoted by (7). The rotated edges E _{ hr } and E _{ vr } of SheppLogan phantom are shown in Figure 3(c) and (f) with the rotation angle of –45°. Based on (5), (7) could be rewritten by the finite differential operators
where D _{ vr } and D _{ hr } are the rotated vertical and horizontal finite differential operators.
MDATV
Based on the multidirection representation, the MDATV is defined as
where E _{ h } and E _{ v } are defined in (5), $\eta $ is the weights adjusting the the energy ratios of the rotated horizontal and vertical edges. Substitute (7) into (8), we get the simplified form of (8)
In this paper we set $\eta $ = 1000 with the same reason in ATV. Therefore, the rotated horizontal edges E _{ hr } are representing the edges parallel to the projection rays. If $\eta $ < 1, then the rotated vertical edges E _{ vr } would represent the edges parallel to the projection rays.
In optimization problem (1), the objective function is some kind of norm function, such as TV, ATV or MDATV. Whatever the objective function is, the optimization requires to minimize this function. One common method for minimization is gradient descent (GD). Thus we need to compute the gradient of MDATV. Substitute (5) into (8) and the MDATV’s gradient is
where
and ϵ is a small positive constant to avoid the singularity. In this paper we set ϵ = 1.0 × 10^{ 8}.
Minimization approaches
GD (Gradient Descent)
The workflow of GD for minimizing regularization in optimization (1) is as following:

a)
Compute the balancing parameter between ART (data fidelity) and GD (regularization)
$$d\phantom{\rule{0.5em}{0ex}}:=\phantom{\rule{0.5em}{0ex}}{{\overrightarrow{f}}^{\left(J\right)}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}{\overrightarrow{f}}^{\left(0\right)}}_{2}$$ 
b)
for k = 1,2,… K
$${\overrightarrow{f}}^{\left(J,0\right)}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}{\overrightarrow{f}}^{\left(J\right)}$$$$\overrightarrow{v}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}\frac{\partial {\overrightarrow{f}}_{\mathit{\text{reg}}}}{\partial {f}_{\mathit{ij}}}$$$$\stackrel{\wedge}{v}}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}\frac{\overrightarrow{v}}{{\left\overrightarrow{v}\right}^{2}$$$${\overrightarrow{f}}^{\left(J,k\right)}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}{\overrightarrow{f}}^{\left(J,k1\right)}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}a\xb7d\xb7\stackrel{\wedge}{v}$$
end for
Where ${\overrightarrow{f}}^{\left(0\right)}$ and ${\overrightarrow{f}}^{\left(J\right)}$ are the estimated image before and after the J th iteration of ART respectively, and $\alpha $ is the step size. Although the variable d is kind of useful for adaptively adjusting the step size of GD, $\alpha $ is determinant in practice.
Note that the norm of $\overrightarrow{v}/{\left\overrightarrow{v}\right}^{2}$ is smaller than the norm of $\overrightarrow{v}/\left\overrightarrow{v}\right$. In GD there usually is a fixed best step size $\alpha $* such that the convergence rate is fastest. This $\alpha $* may be found by the traversal method, which tries many step sizes for the problem and finds the best one within a given accuracy. Obviously, smaller the gradient vector $\stackrel{\wedge}{v}$ is, bigger the best step size $\alpha $* is. The bigger step size is more resistant to the perturbation. For example, the perturbation of 0.1 is negligible for the step size of 50, but this perturbation would have a significant impact for the step size of 0.01. Therefore, we use $\overrightarrow{v}/{\left\overrightarrow{v}\right}^{2}$ instead of $\overrightarrow{v}/\left\overrightarrow{v}\right$ as the gradient vector $\stackrel{\wedge}{v}$ in this paper.
However, in optimization (1), GD is not a stable method for minimizing TV, ATV or MDATV regularizations. Since the best step size of GD changes greatly with the projection parameters and the target images. And the imperfect step size may greatly slow down the convergence rate and degenerate the quality of the estimated image. To overcome this defect of GD, the NESTA (NESTerov’s algorithm) [36] method is proposed to minimize the regularizations in (1).
NESTA
Because the MDATV could be expressed as
where $\mathit{E}{\mathit{\text{'}}}_{\mathit{hr}}\phantom{\rule{0.1em}{0ex}}=\sqrt{\eta}{\mathit{E}}_{\mathit{hr}},\xb7\mathit{E}{\mathit{\text{'}}}_{\mathit{vr}}\phantom{\rule{0.1em}{0ex}}=\phantom{\rule{0.1em}{0ex}}{\mathit{E}}_{\mathit{vr}}$. The form of MDATV norm in (10) is the same as the TV norm in (3), which is almost the similar for ATV and TV. Thus we can first study the NESTA method for TV minimization and then apply the conclusions to ATV and MDATV minimization.
NESTA, based on the Nesterove’s smoothing technique [37], is a fast firstorder method for sparse recovery. First, the nonsmooth TV norm can be rewritten as
where $\overrightarrow{f}\in {\mathcal{Q}}_{p}$, the convex set ${\mathcal{Q}}_{p}$ is referred to as the primal feasible set. u = [u _{1}, u _{2}]^{T} is in the dual feasible set ${\mathcal{Q}}_{d}$ if and only if ${u}_{1}^{2}\left[i,j\right]\phantom{\rule{0.5em}{0ex}}+\phantom{\rule{0.5em}{0ex}}{u}_{2}^{2}\left[i,j\right]\phantom{\rule{0.5em}{0ex}}\le \phantom{\rule{0.5em}{0ex}}1$, and Df _{ i,j } = [D _{ h } f _{ i,j } ,D _{ v } f _{ i,j }]^{T} is the finite difference of the 2D image f _{ i,j }. The superscript T denotes transpose. If u is regarded as a 2D vector, then ${\mathcal{Q}}_{d}$ is a ℓ _{ 2 } norm unit ball.
According to Nesterov’s work, the smoothed TV regularization function is (12)
where $\mu $ should be set sufficiently small such that . Here we set $\mu $ = 0.01. Then is given by(13)
Where D = [D _{ h } ,D _{ v }]^{T}, ${u}_{u}\left(\overrightarrow{f}\right)$ is of the form [u _{ 1 } ,u _{ 2 }]^{T} and for each a ∈ {h, v},
The minimization of the smoothed TV norm by Nesterov’s algorithm is obtained by iteratively estimating three sequences $\left\{{\overrightarrow{f}}_{k}\right\}$, $\left\{{\overrightarrow{d}}_{k}\right\}$, and $\left\{{\overrightarrow{e}}_{k}\right\}$ as follows:
Initialize $\overrightarrow{{f}_{0}}$. for k = 1,2,…,K,

a)
Compute

b)
Compute ${\overrightarrow{d}}_{k}$
(14) 
c)
Compute ${\overrightarrow{e}}_{k}$
(15) 
d)
Update ${\overrightarrow{f}}_{k}$
$${\overrightarrow{f}}_{k+1}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}{\tau}_{k}{\overrightarrow{e}}_{k}\phantom{\rule{0.5em}{0ex}}+\phantom{\rule{0.5em}{0ex}}\left(1{\tau}_{k}\right){\overrightarrow{d}}_{k}$$(16)
end for.
In the above algorithm, α _{ i } = (i + 1)/2, τ _{ k } = 2/(k + 3) is suggested for fast convergence. Since the goal is minimizing TV norm without other constraints, ${\overrightarrow{d}}_{k}$ and ${\overrightarrow{e}}_{k}$ can be computed by letting the gradients of the object functions be zero and then we get
In (15), a suitable smoothing proxfunction is
Then we have
Some other preset parameters are
ART + MDATV scheme
In this paper, the data fidelity constraint is processed by ART iteration, and the MDATV regularization is minimized by NESTA. To use the prior information of the edges in multiple directions, the original ART + TV scheme [11] is not suitable for MDATV regularization. In the ART + TV scheme, firstly, the overall projection data are used to update the estimated image through ART. Then the estimated image is used as the input of the minimization method. The ART and the minimization is operated alternately. Each pair of successive ART and minimization constitutes an iteration of the ART + TV scheme. The workflow of the ART + TV scheme is shown below:
In the above scheme, the stopping criterion is the relative error between the estimated image and the true phantom image being less than 1.0 × 10^{–5}. The relative error is computed as
where ${\overrightarrow{f}}_{r}$ is the estimation of the restored image, ${\overrightarrow{f}}_{0}$ is the original phantom image known as ground truth.  ·  is the Frobenius norm.
However, the MDATV use the prior information of edges in multiple directions by adjusting the angle $\theta $ in (7). In the ART + TV scheme, there is no time to adjust $\theta $. Therefore, the ART + MDATV scheme is proposed.
In the ART + MDATV scheme, projection data are divided into m groups based on the directions of the projection rays. In each group of projections, the projection rays are in the same direction for the parallel projection, and in the approximately same direction for the fanbeam projection. Thus each group of projections corresponds to an angle θ _{ k } (k = 1, 2, …, m), which specifies the location of the Xray source. Then for the k th group of projections, the ART and the minimization of MDATV regularization of θ _{ k } are operated once successively. The estimated image of the k th group of projections is used as the initial estimate of the (k + 1)th group of projections. The process of the m groups of projections constitutes an iteration of the ART + MDATV scheme. The workflow of the ART + MDATV scheme is shown below:
The stopping criterion in the above scheme is the same as that in the ART + TV scheme. Based on the experiments, the iteration times of NESTA for MDATV is set as 10, and the iteration times of NESTA for TV and ATV are both set as 40.
Simulations and experiments
Numerical simulations and real data experiments are conducted to validate the performance of MDATV based approaches. The TV and ATV based approaches are used for comparison. To compare the stableness of GD and NESTA, for each method, they are respectively used for minimizing the regularizations. Therefore, the combination of ART with three different regularizations and two kinds of minimization methods leads to six approaches: ARTTVGD, ARTTVNESTA, ARTATVGD, ARTATVNESTA, ARTMDATVGD, and ARTMDATVNESTA. The simulations for noisy measurements are conducted for comparing the noise robustness of different regularizations. Finally, the real data experiment is used to testify the effectiveness of MDATV in practical applications. Besides, this paper uses the CT module in “Image reconstruction toolbox” [38] to construct the projection matrix and operate the forward and backward projections.
Numerical simulation
Simulation settings
To reduce the Xray dose, reduction of projection rays is a natural option. There are two common styles for reducing the projection rays. They are limitedangle [9] and fewview [10] styles. The schematic diagrams of these two styles are shown in Figure 5. For the same number of projection views, the information of target image in the fewview style is much more than that in the limitedangle style. In this paper, we only describes the simulations of the fewview style.The phantom used in this paper is generated by the function ‘phantom(128)’ in MATLAB^{®}, as shown in Figure 3(a). It has 128 × 128 pixels. The simulated Xray detector has 240 bins for the fanbeam projection. And the size of the detector bin is 2 millimeters. In the fanbeam projection, the distance between the source and center of the detector is 960.45 millimeters, and the distance from the source to the origin is 628.88 millimeters.
According to the ERP (Exact Reconstruction Principle) [6], if the number of FT (Fourier Transform) samples is twice the number of nonzero pixels in the gradient image, then the optimization program (1) can yield a unique solution. The gradient image of SheppLogan phantom is shown in Figure 3(d), where the number of nonzero pixel is 1743. In the digital condition, according to the central slice theory [34], one projection data corresponds to one FT sample. Therefore, we need at least 1743 × 2 = 3486 projection data points. While the fanbeam detector has 240 bins, thus based on ERP we need 3486 ÷ 240 ≈ 15 views of projections. For comparison, we also take some less views of projections, such as 11 and 13 views, in the numerical simulations.
The best step sizes of GD are different when the projection style or the regularization changes. Table 1 lists the best step sizes of GD for the fanbeam projections with different regularizations and projection views. These best step sizes are obtained by the traversal method within the accuracy of 1. Since the perturbation within 1 rarely affects the convergence rate of GD.
To test the robustness of these methods to noise, the Poisson noise is introduced into the projections by the function ‘poisson’ in the “Image reconstruction toolbox” [38]. The incident photon number is set as 1.0 × 10^{5}. Table 2 lists the best step sizes of GD for the noisy reconstructions.
Visualizationbased evaluation
The restored images from 11 views of noisefree and noisy projections are shown in Figure 6, and the restored images from 15 views of noisefree and noisy projections are shown in Figure 7. Both the observations of Figures 6 and 7 indicate that the TV and MDATV regularizations are obviously better than ATV regularization. There are some horizontal artifacts in the restored image of ATV methods, which is caused by the unbalanced regularization, while the projection data are uniformly distributed. It can also be observed that increased projection data prominently improved the image quality.
Profilebased evaluation
To further visualize the difference among various methods, vertical profiles of restored images were drawn across the 64th column, from the 1st row to the 128th row. Figures 8 and 9 shows the vertical profiles of the restored images in Figure 6 and the corresponding profile in the original SheppLogan phantom. And the vertical profiles of the restored images in Figure 7 and their corresponding profile in the original SheppLogan phantom are shown in Figures 10 and 11. The observations from all the profiles also indicate the same conclusion as the visualizationbased evaluation.
Relative error study
The relative error is used in the simulation iterations as stopping criterion. Therefore, it can also be used to show the convergence condition of various methods. The curves of relative errors versus iteration times of reconstruction processes from 11 and 15 views of projections are shown in Figures 12 and 13. The observations indicate that after some iterations the reconstructions would converge to steady conditions. And the relative error after the final iteration can represent the reconstruction quality. For the 11 views condition, MDATV is better than TV, and TV is better than ATV. While for the 15 views condition, the three regularizations are hard to distinguish.
UQI (Universal Quality Index) study
To perform more quantitative analysis of these methods, the UQI [39] is introduced. UQI measures the similarity between the desired image and its ground truth image [23]. UQI value ranges between zero and one. A higher UQI value represents a higher similarity between the testing image and the ground truth image, and vice versa. The ROI (Region Of Interest) used for computing UQI is in the red rectangle as shown in Figure 3(a). The UQI values of restored images by various methods from different views of projections are plotted as several curves shown in Figure 14. The UQI values are in accord with the relative errors. The observations also indicate that the MDATV is better than TV, and TV is better than ATV. When the volume of projection data increase, the distinctions among these regularizations are decreasing.
Real data experiment
To demonstrate the effectiveness of proposed method in real CT application, the real data experiment is conducted. The projection data is from a microCT machine in Tianjin Key Laboratory of Biomedical Detecting Techniques and Instruments. The scanning phantom is an organic glass column with three holes of different diameters, as shown in Figure 15. In this experiment, the left hole on top is empty (filled with air), the right hole on top is filled with corn flour, and the bottom hole is filled with a copper bar. The distance between the source and center of the detector is 960.45 millimeters, and the distance from the source to the origin is 628.88 millimeters.
In the experiment, the detector has 1024 bins, and the physical size of each bin is 0.05 millimeters. Therefore, the raw sinogram of one view has 1024 pixels, and each pixel represents the physical length of 0.05 millimeters. However, the computation load of the iterative reconstruction is very heavy for the high resolution sinogram. Thus the high resolution sinogram is down sampled by the ratio of 4. And the low resolution sinogram of one view has 256 = 1024 ÷ 4 pixels. To maintain the physical length of the projection geometry, each pixel in the low resolution sinogram represents the physical length of 0.2 = 0.05 × 4 millimetres. Actually, in this down sampling, four neighboring pixels in the raw sinogram of one view is averaged as one pixel in the low resolution sinogram of one view.In this experiment, restrained by the microCT machine, the whole projection data contains 120 views of fanbeam projection uniformly distributed across 360°. Its FBP reconstruction is shown in Figure 16. Compare Figure 15 and 16, it is easy to distinguish the air and copper filled holes from the background organic glass, while the corn filled hole is hard to distinguish from the background organic glass. This is because that the attenuation coefficient of corn flour is very similar to that of the organic glass, while the attenuation coefficients of air and copper are very different from that of the organic glass.
For the sparse CT reconstruction experiment, 30 views of fanbeam projections are chose uniformly across the range of 360°. Due to the ATV’s poor performance for the subERP projections (volume of projections is less than that required by ERP), we only compare TV and MDATV based methods. The reconstruction methods used in the experiments are ARTTVGD, ARTMDATVGD, and ARTMDATVNESTA. The step sizes of GD used in the reconstructions are 92 for TV and 10 for MDATV, which are estimated from an approximate synthetic phantom. The iteration times of NESTA for MDATV is 40.In Figure 17, we show the reconstructions after 5 iterations (top row) and 100 iterations (bottom row) with three methods (from left to right): ARTTVGD, ARTMDATVGD, ARTMDATVNESTA. In the restored images, the holes filled with air and copper are distinct, while the hole filled with corn is undistinguishable. To further visualize the difference among various methods, horizontal profiles of restored images were drawn across the 63rd row (for the air and corn filled holes) and 127th row (for the copper filled hole), from the 10th column to the 160th column. The profiles of restored images after 5 and 100 iteration are shown in Figure 18. Since the projection data is not sufficient, the FBP reconstructions abound with artifacts, but the FBP reconstructions clearly depict the boundaries of the air and copper filled holes. And the observations from Figure 18 indicate that MDATVNESTA profiles matches the FBP profiles best.The regions of holes filled with air and copper (in the red rectangles in Figure 17) are zoomed in and shown in Figures 19 and 20 respectively. To increase the contrast of Figures 19 and 20, 5% of data is saturated at low and high intensities of the original image. Figures 19 and 20 indicate that ARTMDATVGD method gave the best results. Since the restored circles by ARTMDATVGD are the most round. The restored circles of ARTMDATVNESTA are more round than that of ARTTVGD. The radial artifacts in the restored images are metal artifacts, which can be removed by some offtheshelf methods.
Discussions and conclusions
The simulations and experiments both indicate that MDATV is a useful and robust regularization for the sparse CT reconstruction when the volume of the projection data is less than the volume required by the ERP. But actually, it is impossible to satisfy the ERP in practical applications. Because the practical target images are not piecewiseconstant, which means the nonzero elements of the gradient image are infinite. Thus the volume of digital projections is always not enough. Therefore, using MDATV to incorporate more a prior information of the target images is valuable and gainful in practical applications, as shown in the real data experiment. And the indistinguishable results among the three regularizations with ERP projections in the simulations is because that the synthetic phantom is piecewiseconstant, which means its ERP samples are limited. Obviously, when the projections’ volume satisfies the ERP, there is no need to incorporate any other prior information.
Since MDATV can represent edges in any direction, the MDATV based methods can use the prior information of edge directions more efficiently than TV or ATV based methods. Therefore, MDATV based methods could enhance the edges in the projection rays directions and suppress the potential artifacts in the no projection rays directions, which leads to the improvement of the reconstructions.
In this paper, NESTA is proposed to minimize the MDATV regularization instead of GD. This is because that NESTA gives almost the same results as GD does for minimization of the regularizations, but NESTA is more stable than GD. For example, Tables 1 and 2 list the best step sizes of GD used in the simulations. The best step size of GD changes much for different projections and noisy conditions. And there is no rule to obtain the best step size of GD in advance. Nevertheless, the few parameters in NESTA almost does not change in the simulations and experiments. Obviously, NESTA is a considerable choice for minimizing the regularizations in the sparse CT reconstruction.
The accumulated computation time of the simulations are listed in Table 3. It is observed that the computation load of NESTA is a little higher than that of GD, but NESTA is more stable than GD. If the comparison includes the time of finding the best step size of GD, then NESTA will be much faster than GD. Due to the special scheme of ART + MDATV, MDATV related methods need more computation time than TV or ATV related methods. This is because that MDATV related methods do the minimization after each group of ART iteration. If the projections are classified into m groups, then MDATV related methods have m  1 times of minimization more than TV or ATV related methods. Furthermore, each minimization has several times of GD or NESTA iterations.
In addition, it is notable that the fanbeam geometry in the numerical simulations and real data experiment are the same, both from the geometry of the microCT machine. In this geometry, the view angle of the fanbeam projection is so small that the projection rays in the same projection view are approximately in the same direction. While in practical applications, the view angle of the fanbeam projection may be bigger so that the projection rays in the same projection view are not approximately in the same direction. In this condition, it needs to rearrange the projection data according to their corresponding projection directions, such that the projection rays having the approximately same direction can be classified into the same group. This will be one of our future research directions.
This paper proposed the MDATV regularization to sufficiently use a prior information of projection directions and image sparsity in the sparse CT reconstruction. Due to the effective use of prior information of projection directions, the restored edge information is greatly enhanced, and thus to improve the reconstructions. The numerical simulations and real data experiments demonstrate the advantage of MDATV over other regularizations. NESTA is proposed as an alternative to GD for minimizing the TVbased regularizations. Because NESTA is more stable and has the almost same performances as GD.
Otherwise, although the MDATV regularization method was only validated for fanbeam projections, it is simple to introduce this regularization approach to other tomography reconstructions. And using prior information of projection directions may further facilitate the development of tomography. Besides, we think the presentation of edges in multiple directions may have more applications in the fields of image enhancement and reconstruction.
Abbreviations
 MDATV:

Multidirection anisotropic total variation
 CT:

Computed tomography
 ATV:

Anisotropic total variation
 TV:

Total variation
 2DIGS:

2Dimensional image gradient space
 ART:

Algebraic reconstruction technique
 NESTA:

NESTerov’s algorithm
 GD:

Gradient descent
 IRIS:

Iterative reconstruction in image space
 ASiR:

Adaptive statistical iterative reconstruction
 GE:

General electric Co
 iDose:

trademark of Philips
 AIDR:

Adaptive iterative dose reduction
 FBP:

Filtered back projection
 ARTTVMIN:

Algebraic reconstruction technique total variation MINimization
 ASDPOCS:

Adaptive steepest descent projection onto convex sets
 PICCS:

Prior image constrained compressed sensing
 PIRPLE:

Prior image registered penalized likelihood estimation
 FCCS:

Feature constrained compressed sensing
 TpV:

Total pvariation
 GPU:

Graphics processing unit
 CG:

Conjugate gradient
 GPBB:

Gradient projection Barzilai Borwein
 ADMM:

Alternating direction method of multipliers
 POCS:

Projection onto convex sets
 ERP:

Exact reconstruction principle
 FT:

Fourier transform
 UQI:

Universal quality index
 ROI:

Region of interest.
References
 1.
Schindera ST, Diedrichsen L, Muller HC, Rusch O, Marin D, Schmidt B, Raupach R, Vock P, SzucsFarkas Z: Iterative reconstruction algorithm for abdominal multidetector CT at different tube voltages: assessment of diagnostic accuracy, image quality, and radiation dose in a phantom study. Radiology 2011,260(2):454–462.
 2.
Martinsen AC, Saether HK, Hol PK, Olsen DR, Skaane P: Iterative reconstruction reduces abdominal CT dose. Eur J Radiol 2012,81(7):1483–1487.
 3.
Deak Z, Grimm JM, Treitl M, Geyer LL, Linsenmaier U, Korner M, Reiser MF, Wirth S: Filtered back projection, adaptive statistical iterative reconstruction, and a modelbased iterative reconstruction in abdominal CT: an experimental clinical study. Radiology 2013,266(1):197–206.
 4.
Wu TH, Hung SC, Sun JY, Lin CJ, Lin CH, Chiu CF, Liu MJ, Teng MM, Guo WY, Chang CY: How far can the radiation dose be lowered in head CT with iterative reconstruction? Analysis of imaging quality and diagnostic accuracy. Eur Radiol 2013,23(9):2612–2621.
 5.
Gervaise A, Osemont B, Lecocq S, Noel A, Micard E, Felblinger J, Blum A: CT image quality improvement using Adaptive Iterative Dose Reduction with widevolume acquisition on 320detector CT. Eur Radiol 2012,22(2):295–301.
 6.
Candès EJ, Romberg J, Tao T: Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory 2006,52(2):489–509.
 7.
Candès EJ, Romberg JK, Tao T: Stable signal recovery from incomplete and inaccurate measurements. Comm Pure Appl Math 2006,59(8):1207–1223.
 8.
Candès EJ, Wakin MB: An introduction to compressive sampling. IEEE Signal Process Mag 2008,25(2):21–30.
 9.
Delaney AH, Bresler Y: Globally convergent edgepreserving regularized reconstruction: an application to limitedangle tomography. IEEE Trans Image Process 1998,7(2):204–221.
 10.
Li M, Yang H, Kudo H: An accurate iterative reconstruction algorithm for sparse objects: application to 3D blood vessel reconstruction from a limited number of projections. Phys Med Biol 2002,47(15):2599–2609.
 11.
Sidky EY, Kao C, Pan X: Accurate image reconstruction from fewviews and limitedangle data in divergentbeam CT. J Xray Sci Technol 2006,14(2):119–139.
 12.
Sidky EY, Pan XC: Image reconstruction in circular conebeam computed tomography by constrained, totalvariation minimization. Phys Med Biol 2008,53(17):4777–4807.
 13.
Chen GH, Tang J, Leng S: Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. Med Phys 2008,35(2):660–663.
 14.
Stayman JW, Otake Y, Prince JL, Khanna AJ, Siewerdsen JH: Modelbased tomographic reconstruction of objects containing known components. IEEE Trans Med Imag 2012,31(10):1837–1848.
 15.
Stayman JW, Dang H, Ding Y, Siewerdsen JH: PIRPLE: a penalizedlikelihood framework for incorporation of prior images in CT reconstruction. Phys Med Biol 2013,58(21):7563.
 16.
Wu D, Li L, Zhang L: Feature constrained compressed sensing CT image reconstruction from incomplete data via robust principal component analysis of the database. Phys Med Biol 2013,58(12):4047–4070.
 17.
Sidky EY, Pan X, Reiser IS, Nishikawa RM, Moore RH, Kopans DB: Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved imagereconstruction algorithms. Med Phys 2009,36(11):4920–4932.
 18.
Yu H, Wang G: SARTtype image reconstruction from a limited number of projections with the sparsity constraint. Int J Biomed Imaging 2010, 2010: 934847.
 19.
Jia X, Dong B, Lou Y, Jiang SB: GPUbased iterative conebeam CT reconstruction using tight frame regularization. Phys Med Biol 2011,56(13):3787–3807.
 20.
Tian Z, Jia X, Yuan K, Pan T, Jiang SB: Lowdose CT reconstruction via edgepreserving total variation regularization. Phys Med Biol 2011,56(18):5949–5967.
 21.
Liu Y, Ma J, Fan Y, Liang Z: Adaptiveweighted total variation minimization for sparse data toward lowdose xray computed tomography image reconstruction. Phys Med Biol 2012,57(23):7923–7956.
 22.
Chang M, Li L, Chen Z, Xiao Y, Zhang L, Wang G: A fewview reweighted sparsity hunting (FRESH) method for CT image reconstruction. J Xray Sci Technol 2013,21(2):161–176.
 23.
Liu Y, Liang Z, Ma J, Lu H, Wang K, Zhang H, Moore W: Total variationstokes strategy for sparseview Xray CT image reconstruction. IEEE T Med Imaging 2014,33(3):749–763.
 24.
Song J, Liu QH, Johnson GA, Badea CT: Sparseness prior based iterative image reconstruction for retrospectively gated cardiac microCT. Med Phys 2007,34(11):4476–4483.
 25.
Jia X, Lou Y, Li R, Song WY, Jiang SB: GPUbased fast cone beam CT reconstruction from undersampled and noisy projection data via total variation. Med Phys 2010,37(4):1757–1760.
 26.
Vandeghinste B, Goossens B, De Beenhouwer J, Pizurica A, Philips W, Vandenberghe S, Staelens S: SplitBregmanbased sparseview CT reconstruction. 11th International meeting on Fully ThreeDimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully 3D 11) 2011, 431–434.
 27.
Park JC, Song B, Kim JS, Park SH, Kim HK, Liu Z, Suh TS, Song WY: Fast compressed sensingbased CBCT reconstruction using BarzilaiBorwein formulation for application to online IGRT. Med Phys 2012,39(3):1207–1217.
 28.
Sidky EY, Jørgensen JH, Pan X: Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle–Pock algorithm. Phys Med Biol 2012,57(10):3065.
 29.
Ramani S, Fessler JA: A splittingbased iterative algorithm for accelerated statistical Xray CT reconstruction. IEEE Trans Med Imaging 2012,31(3):677–688.
 30.
Niu T, Zhu L: Accelerated barrier optimization compressed sensing (ABOCS) reconstruction for conebeam CT: phantom studies. Med Phys 2012,39(7):4588–4598.
 31.
Jin X, Li L, Chen Z, Zhang L, Xing Y: Anisotropic total variation for limitedangle CT reconstruction. Nuclear Science Symposium Conference Record (NSS/MIC), IEEE 2010, 2232–2238.
 32.
Chen Z, Jin X, Li L, Wang G: A limitedangle CT reconstruction method based on anisotropic TV minimization. Phys Med Biol 2013,58(7):2119.
 33.
Quinto ET: Tomographic reconstructions from incomplete datanumerical inversion of the exterior Radon transform. Inverse Probl 1988,4(3):867.
 34.
Bracewell RN: Strip Integration in Radio Astronomy. Aust J Physics 1956,9(2):198–217.
 35.
Goodman JW: Introduction to Fourier optics, 2nd edn: McGrawHill Companies. 1996.
 36.
Becker S, Bobin J, Candès E: NESTA: a fast and accurate firstorder method for sparse recovery. SIAM J Imaging Sci 2011,4(1):1–39.
 37.
Nesterov Y: Smooth minimization of nonsmooth functions. Math Program 2005, 103: 127–152.
 38.
Fessler J: Image reconstruction toolbox. [http://web.eecs.umich.edu/~fessler/code/]
 39.
Wang Z, Bovik AC: A universal image quality index. Signal Processing Letters, IEEE 2002,9(3):81–84.
Acknowledgment
The authors would like to thank Professor Feng Gao and Professor HuiJuan Zhao for providing microCT data.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
HL and XC conceived the study. HL implemented the studies and drafted the manuscript. XC, YW, and DY contributed to discussions and suggestions throughout this work. XC and DY supervised the study. ZZ and QZ designed and performed the experiments. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Li, H., Chen, X., Wang, Y. et al. Sparse CT reconstruction based on multidirection anisotropic total variation (MDATV). BioMed Eng OnLine 13, 92 (2014). https://doi.org/10.1186/1475925X1392
Received:
Accepted:
Published:
Keywords
 Sparse CT
 Iterative reconstruction
 Anisotropic total variation
 Coordinate rotation transform
 Edges in multiple directions