Sparse-view image reconstruction via total absolute curvature combining total variation for X-ray computed tomography

Abstract

Sparse-view imaging is a promising scanning approach which has fast scanning rate and low-radiation dose in X-ray computed tomography (CT). Conventional L1-norm based total variation (TV) has been widely used in image reconstruction since the advent of compressive sensing theory. However, with only the first order information of the image used, the TV often generates dissatisfactory image for some applications. As is widely known, image curvature is among the most important second order features of images and can potentially be applied in image reconstruction for quality improvement. This study incorporates the curvature in the optimization model and proposes a new total absolute curvature (TAC) based reconstruction method. The proposed model contains both total absolute curvature and total variation (TAC-TV), which are intended for better description of the featured complicated image. As for the practical algorithm development, the efficient alternating direction method of multipliers (ADMM) is utilized, which generates a practical and easy-coded algorithm. The TAC-TV iterations mainly contain FFTs, soft-thresholding and projection operations and can be launched on graphics processing unit, which leads to relatively high performance. To evaluate the presented algorithm, both qualitative and quantitative studies were performed using various few view datasets. The results illustrated that the proposed approach yielded better reconstruction quality and satisfied convergence property compared with TV-based methods.

Keywords

Sparse-view image reconstruction variational method total absolute curvature total variation alternating direction method

1 Introduction

X-ray computed tomography (CT) has been widely used in medical diagnosis and industrial nondestructive testing. However, the ionizing radiation with excessive dose can do serious harm to living tissue cells [1]. On the other hand, the data acquisition may be restricted by the condition of scanning situation, such as time or space limitations. Therefore, sparse-view scanning scheme is proposed to tackle this type of reconstruction procedure. However, reconstruction from sparse- view, in fact, is an ill-posed problem which has put forward challenges for the design of reconstruction algorithms. The quality of the CT image construction with the sparse data sampling could be degraded if improper methods are applied. Because of the violation of Tuy-Smith condition [2] in sparse-view data acquisition, iterative techniques are often used to provide better reconstructions than analytical methods [3]. In dealing with ill-posed problems, iterative algorithms based on the variational principle of regularization algorithm are the most common methods. For regularization method, the choice of the appropriate regularization term is the first issue to be considered.

Compressive sensing (CS) [4 –6] has been developed as a powerful tool for exact recovery of signals or images from under-sampled observations. The exact reconstruction is based on the fact that there is a sparse representation of an image in transformed domain, such as image gradient. In this case, total variation (TV) regularization and the L1-norm of frame coefficients is often applied in image reconstructions and denoising applications [7, 8]. However, conventional TV has several unfavorable features, such as stair-case effect, that is, the recovery image has piecewise-constant blocks even if the original image is smooth [9]. To eliminate stair-case effect, the straightforward way is involving higher order image derivative into regularizer to provide better edge preservation while reducing the staircase effect in smooth regions of image [10]. In previous work, higher order regularization models which constitute valid extensions of TV are effectively introduced for more image processing tasks [11, 12]. However, higher order regularization models based on second order derivatives are difficult to solve the numerical problems with fourth-order partial differential equations (PDE). Besides, the results generated by some fourth-order problems may still look patchy in smooth regions [10].

It is also highly desirable to provide such regularization term with strong geometric foundations in such a way that mathematical tools from the fields of geometry can be properly analyzed. Recently, there has been some improvement in this field with the development of geometric entities such as the total curvature (TC) approach for surface fairing [13], the mean curvature (MC) method for imaging denoising [14] and the work of Carlos and Chen using Gaussian curvature (GC) for denoising [15]. Applying the L1-norm, TV regularization term merely utilizes the sparsity of gradient (first-order information). Beyond the gradient information, there are many other geometric attributes, such as curvature, which can serve as characteristic description of the featured complicated images. In this paper, the curvature information with the L1-norm is adapted into regularization term to penalize sparsity more effectively and a new hybrid model containing TAC and TV of the image surface is proposed.

For variational model containing the second derivative, it is a nontrivial issue to construct efficient numerical algorithms to deal with the corresponding fourth order Euler-Lagrangian equation. In the recent years, there have been some improvements in the development of algorithms such as the nonlinear multi-grid method of Chumchob [16] and the two-step (TS) algorithm of Tasdizen [17], the convexity splitting method [18] and the fixed-point homotopy method of Yang [19].

There is also the Augmented Lagrangian (AL) method which has already been successfully utilized in the variational models in image processing [20, 21]. Tai and Hahn 20] proposed a fast and efficient algorithm to solve the minimization problem related to Euler’s Elastica (EE) energy. This method decomposed the original minimization problem into a few simple ones, some of which can be solved by the FFT, while the other ones can be quickly dealt with due to existing explicit solutions. An alternating direction total variation minimization (TV-ADM) algorithm [22] for sparse-view CT reconstruction has been proposed, using alternating direction method (ADM) based on ALM. Wang [23] proposed an inexact-TV-ADM algorithm which efficiently computed the pseudo inverse and hence reduced the computing time with little loss of accuracy. Inspire by [20, 21], in this paper, we propose a fast and efficient algorithm to solve the minimization problem related to TAC-TV energy.

The rest of this paper is organized as follows. In section 2, the proposed model is introduced in detail. We review the variational geometry and propose a total absolute curvature combining total variation model. A detailed procedure is discussed for solving this model utilizing ALM, and numerical implementation is also developed. To evaluate performance of the proposed method, we also propose the parameter selection strategy, stopping criterion based on KKT conditions, and texture distance of ROI regions. In section 3, the numerical results are demonstrated to test the proposed method. We list numerical experiments for our model, and compare our results with those obtained by the alternating direction total variation minimization (TV-ADM). Finally in section 4, we present our discussion and conclusion.

2 Method

2.1 Some backgrounds for variational methods and curvature-based regularization

In computed tomography image reconstruction, conventionally, satisfying images are obtained from a large number of projection views [2]. It is of practical value to develop imaging techniques capable of yielding accurate reconstructions from sparse projections (referred to as sparse-view data). However, image reconstructions from few-view projections are actually ill-posed inverse problems, which generally are quite not easy to find the satisfying solutions, especially due to the lack of sufficient measurements and the presence of noise in the observation.

To this end, variational models are often taken to address the above mentioned concerns in both image processing and CT image reconstructions. Such models are built on the minimization of a well-chosen energy function, typically as $u^{*} = arg min_{u} {H (u) + E (u)},$ where H (u) is a fidelity term which ensures that the estimated value is consistent with the observation data. It is noted that the quadratic term based on L2-norm is often applied as the fidelity term. The second term E (u) is a regularization term which is utilized to describe some critical properties of the intended solution. For better reconstruction, a suitable regularization term E (u) is of crucial importance and this paper is trying to make some efforts to establish efficient forms.

One of the most popular variational model is proposed by Rudin, Osher, and Fatemi (ROF) [7], known as the total variation (TV) methods which are first proposed to recover images contaminated by noise. Although its mathematical expression is simple and intuitive, yet it can be explained from the view of differential geometry. Let u (x) be a parametric curve L defined on x ∈ Ω. A first natural geometric measure is the length [24]: $| L | : = E_{1} = \int_{L} ds = \int_{Ω} | \nabla u | dx,$ where the length |L| can be interpreted as the total variation of u, and ∇u is the gradient of this curve. The definition of TV can be extended to multi-dimensional signals straightforwardly, such as images and videos. The remarkable feature of the ROF model lies in its effectiveness in preserving object edges while removing noise. The TV-based method has been widely employed to accomplish other image processing tasks such as denoising, deblurring and segmentation.

Moreover, since the advent of compressed sensing, TV-based reconstructions have become the mainstream of the research focus for CT image reconstruction. Conventional TV-based methods apply the first order information of the images in the form of L1-norm. Although the ROF model has been studied and utilized popularly, but it still has several unfavorable features such as stair-case effect, corner smearing and loss of image contrast. It is known that the ROF model contains only the first order information. Therefore, in order to remedy these unfavorable properties of this model, one can consider incorporating second order derivative information into energy term E (u).

Curves and surfaces are basic geometric elements in image and vision analysis, and computer graphics. In vision psychology, it has been observed that human subjects can easily observe and process convex or concave features of shapes and boundaries [24], which can be conveniently described by the curvature. That is to say, curvature plays an important role in biological vision networks. In two dimensional situations, for instance, a general second order measure based on curvature measurement for a smooth curve is given by $E = \int_{Ω} ω (| κ |) {| κ |}^{p} ds,$ where |κ| is the absolute value of curvature, ω is nonnegative and increasing with respect to |κ|, and constant p ≥ 1. A specific and famous form is Euler’s elastic (EE) energy: $E_{EE} = \int_{Ω} (a + b κ^{2}) ds,$ where a and b are positive constants. The EE model was first introduced in the field of computer vision by Mumford [25] for measuring the quality of interpolating curves in disocclusion. To some extent, although the EE model is very effective in many areas of image processing, it cannot take sufficient advantage of the sparsity for low-dose CT image reconstruction from sparse-view data. For the smoothly changing voxel values in the desired image, such as gradual change regions, their gradient magnitude image (GMI) is not sparse, but the curvature of the image may contain some sparse information. Thus, by constructing an appropriate regularization function to describe the sparsity of curvature, it may be possible to provide a feasible priori for sparse problem. In order to penalize sparsity more effectively, one may weights the angular information, then the total absolute curvature (TAC) $\int_{Ω} | κ | ds = \int_{Ω} | \frac{d θ}{ds} | ds,$ is exactly the total variation of the heading direction θ. The TAC applies L1-norm as the regularization term, which has taken the sparsity of the curvature (second order information) into consideration. In this paper, we take advantage of the sparsity of curvature and the length simultaneously. Then the total absolute curvature combing total variation energy can be expressed as $E^{*} = \int_{Ω} (a + b | κ |) ds,$ where a and b are some positive constant weights, the ratio a/b indicates the relative importance of the total length versus total squared curvature.

For a function u (x), the curvature of the level curve l_c : = {x|u (x) = c} is $| κ | = | \nabla \cdot \frac{\nabla u}{| \nabla u |} |$ . Thus, the TAC-TV energy for the level curve l_c is given by $E_{l_{c}} = \int_{l_{c}} [a + b | \nabla \cdot \frac{\nabla u}{| \nabla u |} |] ds .$

Using the co-area formula to integrate all the curves [24], the TAC-TV energy is as follows: $E^{*} (u) = \int_{Ω} [a + b | \nabla \cdot \frac{\nabla u}{| \nabla u |} |] | \nabla u | dx,$ where ∇ denotes the differential operator, and dx = dx₁dx₂.

2.2 TAC-TV minimization for CT Image reconstructions

Under the discrete-to-discrete model, the computed tomography (CT) imaging model can be expressed using the following discrete linear system: $Au = b,$ where $u \in ℝ^{N}$ is the image vector comprised of the voxel coefficients, $A \in ℝ^{M \times N}$ is the system matrix generated by projection operation, and $b \in ℝ^{M}$ represents the observed data vector containing the measured projection data. Sparse-view projection data are not sufficient for exact reconstruction by conventional methods, which means that the equation is ill-posed and the unique solution is difficult to obtain. Conventional TV minimization is useful and practical in computational aspects but sacrifices reconstruction capability for extremely under-sampled data.

As discussed in the above context, it can be pointed out that TAC-TV energy is a more suitable model for featured complicated CT image reconstruction than TV term. Therefore, we proposed a hybrid model containing both TV and the TAC for obtaining the intended CT image, termed as TAC-TV optimization which can be expressed as: $u^{*} = arg min_{u} (a + b | \nabla \cdot \frac{\nabla u}{| \nabla u |} |) | \nabla u |, s . t . Au = b,$ (1) where | ∇ u| = ∑_i ∥ D_iu ∥ ₁, D_i denotes the finite differential operator along direction i. Specifically, D₁ and D₂ denote the horizontal and vertical differential operator respectively in two dimensional form. In this paper, in convenience, it should be emphasize here that when a scalar plus/minus a vector, it means all the entities in this vector should plus/minus this scalar, returning a vector and vice versa. For instance, 〈a, u〉 is equal to 〈a1, u〉 where 1 denotes the row vector with all entries equal to one.

2.3 The Augmented Lagrangian method to the TAC-TV model

As is well known, the augmented Lagrangian method (ALM) has been successfully applied to the constraint optimization problem. Furthermore, we apply the ALM to tackle with the above optimization. Specifically speaking, we first convert the optimization problem (1) as a multi-constraint optimization problem by introducing two new variables p = ∇ u and $n = \frac{\nabla u}{| \nabla u |}$ : $\begin{matrix} u^{*} = arg min_{u} (a + b | \nabla \cdot n |) | p |, \\ s . t . Au = b, p = \nabla u, n = \frac{\nabla u}{| \nabla u |} . \end{matrix}$ (2)

An unconstrained form of the above optimization problem, using the ALM, can be expressed as follows: $\begin{matrix} L (u, p, n; λ_{0}, λ_{1}, λ_{2}) = (a + b | \nabla \cdot n |) | p | + \frac{r_{0}}{2} {∥ Au - b ∥}^{2} + < λ_{0}, Au - b > \\ + \frac{r_{1}}{2} {∥ | p | n - p ∥}^{2} + < λ_{1}, | p | n - p > \\ + \frac{r_{2}}{2} {∥ p - \nabla u ∥}^{2} + < λ_{2}, p - \nabla u >, \end{matrix}$ (3) where 〈, 〉 denotes inner product and · denotes dot product.

However, inspired by previous studies [20, 21], instead of dealing with optimization problem as (1.3), we here consider a more relaxed and easily handled problem. Beforehand, the following useful results are given in [21], and for the self-completeness we include them here:

Lemma 1. Given two vectors n ≠ 0, |p| ≠ 0. If they satisfy $| n | \leq 1, | p | = n \cdot p,$

Then we have $n = \frac{p}{| p |} .$

As a consequence of Lemma 1, it is easy to see that minimization of (1) is equivalent to the following constrained minimization problem by introducing variables p, m, n and v: $\begin{matrix} u^{*} = arg min_{u} (a + b | v |) | p |, \\ s . t . Au = b, p = \nabla u, | p | = m \cdot p, | m | \leq 1, n = m, v = \nabla \cdot n, \end{matrix}$ (4) where a, b > 0, $u, v \in ℝ^{N},$ $p, m, n \in ℝ^{N \times 2},$ $b \in ℝ^{M},$ $A \in ℝ^{M \times N}$ .

The corresponding augmented Lagrangian function can be formulated as $\begin{matrix} L (u, p, m, n, v; λ_{0} {, λ}_{1} {, λ}_{2} {, λ}_{3} {, λ}_{4}) = (a + b | v |) | p | \\ + \frac{r_{0}}{2} {∥ Au - b ∥}^{2} + < λ_{0}, Au - b > \\ + \frac{r_{1}}{2} {∥ p - \nabla u ∥}^{2} + < λ_{1}, p - \nabla u > \\ + r_{2} (| p | - m \cdot p) + λ_{2} (| p | - m \cdot p) + δ_{R} (m) \\ + \frac{r_{3}}{2} {∥ n - m ∥}^{2} + < λ_{3}, n - m > \\ + \frac{r_{4}}{2} {∥ v - \nabla \cdot n ∥}^{2} + < λ_{4}, v - \nabla \cdot n >, \end{matrix}$ (5) where $R = {m \in L^{2} (Ω) : | m | \leq 1 a . e . in Ω}$ and the indicator function on $R$ and can be expressed as follows: $δ_{R} (m) = {\begin{matrix} 0, & | m | \leq 1; \\ + \infty, & otherwise . \end{matrix}$ (6)

The introduction of the variable m is to relax the variable n, and the variable m is required to lie in the set $R$ so that the term |p| - m · p is always non-negative. The non-negativity is the reason why we use |p| - m · p as a penalization instead of L2 penalization. The application of ALM leads the following algorithm framework, listed in Table 1, for the TAC-TV minimization for CT image reconstructions.

Table 1

Augmented Lagrangian method for the TAC-TV model

1. Initialization: u⁰, p⁰, m⁰, n⁰, v⁰, and

λ_{0}^{0}, λ_{1}^{0}, λ_{2}^{0}, λ_{3}^{0}, λ_{4}^{0} .

2. For k ≥ 1, compute an approximate minimizer (u^k, p^k, m^k, n^k, v^k) of the augmented

Lagrangian function with the fixed Lagrangian multiplier

λ_{0}^{k - 1}, λ_{1}^{k - 1}, λ_{2}^{k - 1}, λ_{3}^{k - 1}, λ_{4}^{k - 1}

(u^{k} {, p}^{k} {, m}^{k} {, n}^{k} {, v}^{k}) \approx L ({u, m, p, n, v; λ}_{0}^{k - 1} {, λ}_{1}^{k - 1}, λ_{2}^{k - 1} {, λ}_{3}^{k - 1} {, λ}_{4}^{k - 1})

3. Update the Lagrangian multipliers

λ_{0}^{k} = λ_{0}^{k - 1} + r_{0} ({Au}^{k} - b),

λ_{1}^{k} = λ_{1}^{k - 1} + r_{1} (p^{k} - \nabla u^{k}),

λ_{2}^{k} = λ_{2}^{k - 1} + r_{2} (| p^{k} | - m^{k} \cdot p^{k}),

λ_{3}^{k} = λ_{3}^{k - 1} + r_{3} (n^{k} - m^{k}),

λ_{4}^{k} = λ_{4}^{k - 1} + r_{4} (v^{k} - \nabla \cdot n^{k}),

4. Measure the relative residuals until a stopping criterion is satisfied

2.4 The reconstruction algorithm by alternating direction method

The augmented Lagrangian function (5) can be split into five sub-problems with respect to u, p, n, m and v, respectively. Therefore, the framework to design the algorithm is based on ADM which is originated to solve a class of problems with multi-variables. ADM framework is very suitable for solving the proposed model in that the proposed model has five separate variables. Moreover, the convergence analysis and conclusions [26] guarantee that this kind of methodology is often of relatively high efficiency. By applying ADM, we need to consider the following sub-problems and obtain the solutions for each of them: $ɛ_{1} (u) = \frac{r_{0}}{2} {∥ Au - b ∥}^{2} + < λ_{0}, Au - b > + \frac{r_{1}}{2} {∥ p - \nabla u ∥}^{2} + < λ_{1}, p - \nabla u >,$ (7) $ɛ_{2} (p) = (a + b | v |) | p | + \frac{r_{1}}{2} {∥ p - \nabla u ∥}^{2} + λ_{1} (p - \nabla u) + (r_{2} + λ_{2}) (| p | - m \cdot p),$ (8) $ɛ_{3} (m) = - (r_{2} + λ_{2}) m \cdot p + \frac{r_{3}}{2} {∥ n - m ∥}^{2} + < λ_{3}, n - m > + δ_{R} (m),$ (9) $ɛ_{4} (n) = \frac{r_{3}}{2} {∥ n - m ∥}^{2} + < λ_{3}, n - m > + \frac{r_{4}}{2} {∥ v - \nabla \cdot n ∥}^{2} + < λ_{4}, v - \nabla \cdot n >,$ (10) $ɛ_{5} (v) = b | v | | p | + \frac{r_{4}}{2} {∥ v - \nabla \cdot n ∥}^{2} + < λ_{4}, v - \nabla \cdot n > .$ (11)

In the following subsections, we will discuss the minimizers of all the above sub-problems and the update formula of the multipliers. In kth iteration for example, for each variable in (4), we fix all the other variables and seek the solution for them to update the corresponding variable. Once all the variables are updated, the multipliers will also be updated. All the process will be repeated until the variables are converged.

A) sub-problem for u

The u sub-problem can be expressed as $\begin{matrix} ɛ_{1} (u) = \frac{r_{0}}{2} {∥ Au - b ∥}^{2} + < λ_{0}, Au - b > + \frac{r_{1}}{2} {∥ p - \nabla u ∥}^{2} + < λ_{1}, p - \nabla u >, \\ = \frac{r_{0}}{2} {∥ A u - b + \frac{λ_{0}}{r_{0}} ∥}^{2} + \frac{r_{1}}{2} {∥ p - \nabla u + \frac{λ_{1}}{r_{1}} ∥}^{2} + C_{1}, \end{matrix}$ (12) where C₁ is a constant independent of the variable u. In order to avoid computing the inverse of the matrix containing A. ${∥ Au - b + \frac{λ_{0}}{r_{0}} ∥}^{2}$ can be linearized at current point u^k: ${∥ Au - b + \frac{λ_{0}}{r_{0}} ∥}^{2} \approx {∥ {Au}^{k} - b + \frac{λ_{0}}{r_{0}} ∥}^{2} + 2 g^{k} (u - u^{k}) + \frac{1}{τ} {∥ u - u^{k} ∥}^{2},$ (13) where τ > 0 and $g^{k} = A^{T} ({Au}^{k} - b + \frac{λ_{0}}{r_{0}})$ is the gradient of ${∥ Au - b + \frac{λ_{0}}{r_{0}} ∥}^{2}$ at u^k. Then (13) transformed into the following form: $ɛ_{1} (u) = r_{0} g^{k} (u - u^{k}) + \frac{r_{0}}{2 τ} {∥ u - u^{k} ∥}^{2} + \frac{r_{1}}{2} {∥ p - \nabla u + \frac{λ_{1}}{r_{1}} ∥}^{2} + \tilde{C_{1}},$ (14)

By setting the derivative to 0, we get $0 = r_{0} g^{k} + \frac{r_{0}}{τ} (u - u^{k}) + r_{1} \nabla^{T} (\nabla u - p - \frac{λ_{1}}{r_{1}}),$ (15) thus $(r_{1} \nabla^{T} \nabla + \frac{r_{0}}{τ} I) u = c_{k},$ (16) where $c_{k} = r_{1} \nabla^{T} p + \nabla^{T} λ_{1} + \frac{r_{0}}{τ} u^{k} - r_{0} g^{k}$ .

Note that, under the periodic boundary conditions for ∇, the coefficient $r_{1} \nabla^{T} \nabla + \frac{r_{0}}{τ} I$ is a constant and block-circulant matrix that can easily be diagonalized by fast Fourier transform (FFT). Let $J = diag [𝔽^{- 1} (r_{1} \nabla^{T} \nabla + \frac{r_{0}}{τ} I) 𝔽]$ , where $𝔽$ stands for the Fourier transform matrix (implemented by FFT) and diag returns a vector constructed by the principal diagonal entities of M. The solution of (16) can be expressed by $u^{k + 1} = 𝔽 [𝔽^{- 1} (c_{k}) / J],$ where the division is a component-wise operator and J is a constant vector that can be computed beforehand to save time.

To obtain the minimization of ɛ₂ (p), ɛ₃ (m), ɛ₄ (n) and ɛ₅ (v). The following lemmas [21] are first introduced.

Lemma 2. Let x be a given scalar or vector, α, β be two parameters with β > 0, then the minimization of the following function $min_{x} α | x | + \frac{β}{2} {∥ x - x_{0} ∥}^{2},$ reads $x = max {| x_{0} | - \frac{α}{β}, 0} sgn (| x_{0} |) .$

Lemma 3. (Lions and Mercier [27]) Let m₀ be a given vector, then the minimization problem $min_{m} {∥ m - m_{0} ∥}^{2} + δ_{R} (m),$ has an explicit solution as $m = {Proj}_{R} (m_{0}) = {\begin{matrix} m_{0}, & | m_{0} | \leq 1; \\ \frac{m_{0}}{| m_{0} |}, & otherwise . \end{matrix}$

B) sub-problem for p

Note that $\begin{matrix} ɛ_{2} (p) = (a + b | v | + r_{2} + λ_{2}) | p | + \frac{r_{1}}{2} {∥ p - \nabla u ∥}^{2} + < λ_{1} + λ_{2} m + r_{2} m, p - \nabla u > \\ = 〈 a + b | v | + r_{2} + λ_{2}, | p | 〉 + \frac{r_{1}}{2} {∥ p - \nabla u + \frac{λ_{1} + λ_{2} m + r_{2} m}{r_{1}} ∥}^{2} + C_{2}, \end{matrix}$ (17) where C₂ is a constant independent of the variable p. We get the minimization of function ɛ₂ (p) by apply the lemma 2 $p = max {| \nabla u - \frac{λ_{1} + λ_{2} m + r_{2} m}{r_{1}} | - \frac{a + b | v | + r_{2} + λ_{2}}{r_{1}}, 0} sgn (\nabla u - \frac{λ_{1} + λ_{2} m + r_{2} m}{r_{1}}) .$ (18)

C) sub-problem for m

The m sub-problem can be express as follows $\begin{matrix} ɛ_{3} (m) = - (r_{2} + λ_{2}) p \cdot m + \frac{r_{3}}{2} {∥ n - m ∥}^{2} - < λ_{3}, n - m > + δ_{R} (m), \\ = \frac{r_{3}}{2} {∥ n + \frac{λ_{3} + (r_{2} + λ_{2}) p}{r_{3}} - m ∥}^{2} + δ_{R} (m) + C_{3}, \end{matrix}$ (19) where C₃ is a constant independent of the variable m. The minimizer of ɛ₃ (m) is provided by the lemma 3 $m = {Proj}_{R} (n + \frac{λ_{3} + (r_{2} + λ_{2}) p}{r_{3}}) .$ (20)

D) sub-problem for n

As for the sub-problem ɛ₄ (n), we expand the variables of n = (n₁, n₂), m = (m₁, m₂), λ₃ = (λ₃₁, λ₃₂) resulting $\begin{matrix} ɛ_{4} (n) = \frac{r_{4}}{2} {∥ \frac{λ_{4}}{r_{4}} + v - \nabla \cdot n ∥}^{2} + \frac{r_{3}}{2} {∥ n - m ∥}^{2} + < λ_{3}, n >, \\ = \frac{r_{4}}{2} {∥ v - D_{1} n_{1} - D_{2} n_{2} + \frac{λ_{4}}{r_{4}} ∥}^{2} \\ + \frac{r_{3}}{2} {∥ n_{1} - m_{1} ∥}^{2} + \frac{r_{3}}{2} {∥ n_{2} - m_{2} ∥}^{2} + < λ_{31}, n_{1} > + < λ_{32}, n_{2} > + C_{4}, \end{matrix}$ (21)

where C₄ is a constant independent of the variable n. We need to get the associated Euler-Lagrange equation. Standard procedures lead to the following equations of n. $\begin{matrix} M_{11} n_{1} + M_{12} n_{2} = N_{1}, \\ M_{21} n_{1} + M_{22} n_{2} = N_{2}, \end{matrix}$ (22) where $M_{11} = r_{4} D_{1}^{T} D_{1} + r_{3},$ $M_{12} = r_{4} D_{1}^{T} D_{2}$ , $M_{21} = r_{4} D_{2}^{T} D_{1},$ $M_{22} = r_{4} D_{2}^{T} D_{2} + r_{3},$ $N_{1} = D_{1}^{T} (r_{4} v + λ_{4}) + r_{3} m_{1} - λ_{31}$ , $N_{2} = D_{2}^{T} (r_{4} v + λ_{4}) + r_{3} m_{2} - λ_{32}$ . Note that, by the properties of difference operator for D_i, the binary linear equations can be easily solved by the 2D-FFT, and here we give the solutions expressed as $(\begin{matrix} λ_{11} & λ_{12} \\ λ_{21} & λ_{22} \end{matrix}) (\begin{matrix} {\tilde{n}}_{1} \\ {\tilde{n}}_{2} \end{matrix}) = (\begin{matrix} {\tilde{N}}_{1} \\ {\tilde{N}}_{2} \end{matrix}),$ (23) $\begin{matrix} n_{1} = 𝔽^{- 1} {\frac{λ_{22} {\tilde{N}}_{1} - λ_{12} {\tilde{N}}_{2}}{λ_{11} λ_{22} - λ_{12} λ_{21}}}, \\ n_{2} = 𝔽^{- 1} {\frac{λ_{11} {\tilde{N}}_{2} - λ_{21} {\tilde{N}}_{1}}{λ_{11} λ_{22} - λ_{12} λ_{21}}} . \end{matrix}$ (24)

E) sub-problem for v

The v sub-problem can be express as follows $ɛ_{5} (v) = b | v | | p | + \frac{r_{4}}{2} {∥ v - \nabla \cdot n + \frac{λ_{4}}{r_{4}} ∥}^{2} + C_{5},$ (25) where C₅ is a constant independent of the variable v. We get the minimization of sub-problem ɛ₅ (v) by apply the lemma 2 $v = max {| \nabla \cdot n - \frac{λ_{4}}{r_{4}} | - \frac{b | p |}{r_{4}}, 0} \cdot sgn (\nabla \cdot n - \frac{λ_{4}}{r_{4}}) .$ (26)

The relationship between the proposed TAC-TV minimization and the conventional TV one is clear: the proposed TAC-TV model includes the TV model, TAC-TV turns to TV reconstruction when b = 0. In this paper, a more relaxed technique is used and the augmented Lagrangian function (5) is transformed into five sub-problems which are iterated in ADM scheme. The general algorithm is summarized in Table 2. In this paper, the proposed method is termed as TAC-TV algorithm.

Table 2

ADM to solve the sub-problems

Initialization: $u_{}^{0} = u^{k - 1}, p_{}^{0} = p^{k - 1}, m_{}^{0} = m^{k - 1}, n_{}^{0} = n^{k - 1},$ and $v_{*}^{0} = v^{k - 1} .$ For l = 1, …, L and fixed Lagrangian
multiplier $λ_{0} = λ_{0}^{k - 1}, λ_{1} = λ_{1}^{k - 1}, λ_{2} = λ_{2}^{k - 1},$ $λ_{3} = λ_{3}^{k - 1},$ and $λ_{4} = λ_{4}^{k - 1},$ compute an approximate minimizer $(u_{}^{l}, p_{}^{l}, m_{*}^{l},$
$n_{}^{l}, v_{}^{l})$ of the augmented Lagrangian function with the fixed Lagrangian multipliers $λ_{0} = λ_{0}^{k - 1}, λ_{1} = λ_{1}^{k - 1}, λ_{2} = λ_{2}^{k - 1},$
$λ_{3} = λ_{3}^{k - 1}, λ_{4}$ $= λ_{4}^{k - 1}$ .
1. Update u(Linearized techniques and 2D FFT)
$u^{k} = 𝔽 [𝔽^{- 1} (c_{k}) / J],$
2. Update p(Shrinkage operator)
$p^{k} = max {\| \nabla u^{k} - \frac{λ_{1} + λ_{2} + r_{2} m^{k - 1}}{} r_{1} \| - \frac{a + b \| v^{k - 1} \| + r_{2} + λ_{3}}{r_{1}}, 0} sgn (\nabla u^{k} - \frac{λ_{1} + λ_{2} + r_{2} m^{k - 1}}{r_{1}}) .$
3. Update m(Projection onto unitary ball)
$m^{k} = {Proj}_{R} (n^{k} + \frac{λ_{3} + (r_{2} + λ_{2}) p^{k}}{r_{3}}) .$
4. Update n(by 2D FFT)
$n_{1}^{k} = 𝔽^{- 1} {\frac{λ_{22} {\tilde{N}}_{1} - λ_{12} {\tilde{N}}_{2}}{λ_{11} λ_{22} - λ_{12} λ_{21}}}, n_{2}^{k} = 𝔽^{- 1} {\frac{λ_{11} {\tilde{N}}_{2} - λ_{21} {\tilde{N}}_{1}}{λ_{11} λ_{22} - λ_{12} λ_{21}}},$
5. Update v(Shrinkage operator)
$v^{k} = max {\| \nabla \cdot n^{k} - \frac{λ_{4}}{r_{4}} \| - \frac{b \| p^{k} \|}{r_{4}}, 0} \cdot sgn (\nabla \cdot n^{k} - \frac{λ_{4}}{r_{4}}) .$
6. Update the multipliers
$\begin{matrix} λ_{0}^{k} = λ_{0}^{k - 1} + r_{0} ({Au}^{k} - b), \\ λ_{1}^{k} = λ_{1}^{k - 1} + r_{1} (p^{k} - \nabla u^{k}), \\ λ_{2}^{k} = λ_{2}^{k - 1} + r_{2} (\| p^{k} \| - m^{k} \cdot p^{k}), \\ λ_{3}^{k} = λ_{3}^{k - 1} + r_{3} (n^{k} - m^{k}), \\ λ_{4}^{k} = λ_{4}^{k - 1} + r_{4} (v^{k} - \nabla \cdot n^{k}), \end{matrix}$
The desired image was obtained by performing the above five procedures in an alternating fashion until the stopping
criterion is satisfied
$(u^{k}, p^{k}, m^{k}, n^{k}, v^{k}) = (u_{}^{l}, p_{}^{l}, m_{}^{l}, n_{}^{l}, v_{*}^{l})$

2.5 Parameter selection and stopping criterion

The augmented Lagrangian penalty parameters r₀, r₁, r₂, r₃ and r₄ are used to balance the data fidelity and four regularization terms. Actually, they are often empirically selected by visual inspection. To get the optimal performance, the values of them should be set in accordance with both the noise level and the sparsity level of the underlying image. Generally speaking, the higher the noise level is, the smaller r₀ should be. In order to guarantee the reconstruction image u is accurate, the value of parameter r₀ could be given larger than that of r₁, r₂, r₃ and r₄, and r₀ ∈ [10², 2 × 10³] is suitable for most applications. To understand the role of the parameter r₁, we go to the equation (16). The parameter r₁ control the amount of diffusion of image u: the larger r₁ is, the more diffusion u has. Generally, the value range of r₁ is [10^-1, 2 × 10²]. It should be noticed that the vector m is introduced to normalize the vector p. So parameter r₂ should be smaller than r₁, and r₁ ∈ [10^-3, 10^-1]. The parameter r₃ controls the amount of closeness between n and m by equation (21), and r₃ ∈ [1, 10³] is reasonable. The parameter r₄ is in [1, 10³] and controls the amount of curvature of image: the larger r₄ is, the more texture u has. While using these as a guideline, we tried to keep the same parameters whenever possible. In particular, our results for the images corrupted by Poisson noise are obtained using the same parameters.

It is important to identify a reasonable and effective stopping criterion for the iteration listed in Table 2 and adapting strategy for the regularization parameter selection in ADMM algorithm. Based on the procedure used in the basic ADMM algorithm [28], we will carry out some theoretical analysis based on the Karush-Kuhn-Tucker (KKT) conditions.

The KKT conditions are satisfied by optimal solutions of constrained optimization problems. Based on (1), an optimal solution $(u^{*}, λ_{0}^{*})$ , commonly called as KKT points, must satisfy the primal feasibility condition: $- A^{T} λ_{0}^{*} \in \partial f (u^{*}),$ (27) here $f (u) : = (a + b | \nabla \cdot \frac{\nabla u}{| \nabla u |} |) | \nabla u |$ , and ∂f denotes the sub-gradient of f. And the dual feasibility condition: ${Au}^{*} - b = 0 .$ (28)

Based on (28) and (29), the residuals of the iteration sequence (u^k, p^k, m^k, n^k, v^k) generated by the algorithm for the feasibility conditions can be used to measure the proximity of the iteration to an optimal solution of the original minimization problem. It is easy to show that the optimality condition of (5): $r_{0} g^{k} + \frac{r_{0}}{τ} (u^{k + 1} - u^{k}) + r_{1} \nabla^{T} (\nabla u^{k + 1} - p^{k} - \frac{λ_{1}^{k}}{r_{1}}) \in \partial f (u^{k + 1}),$ (29) where $g^{k} = A^{T} ({Au}^{k} - b + \frac{λ_{0}^{k}}{r_{0}})$ . The optimality condition leads to that the following quantity which can be viewed as residual for the feasibility conditions in (5). $R_{u}^{k + 1} = r_{0} (A^{T} A - \frac{1}{τ} I) (u^{k + 1} - u^{k}) - r_{1} \nabla^{T} (\nabla u^{k + 1} - p^{k} - \frac{λ_{1}^{k}}{r_{1}}),$ (30) where $R_{u}^{k + 1}$ denotes the k+1th residual of variable u.

In the same manner the optimality condition for the dual feasibility conditions is: $R_{λ_{0}}^{k + 1} = A u^{k + 1} - b,$ (31)

$R_{λ_{0}}^{k + 1}$ denotes the k+1th residual of variable λ₀.

Based on relations (30) and (31), the iteration sequence in the ADMM should be terminated when the residuals for the feasibility conditions (at least one i.e. the smaller of the two) should be small in an predetermined value.

The residuals also offer a strategy to adjust the AL penalty parameters: r₀ and r₁. Considering the AL function, it can be observed that large value of r₀ will induce a small primal residual of λ₀, but will violate the dual residual of u by its definition. Conversely, the definition of residual of u suggests that small value of r₀ will reduce the dual residual, which will also reduce the regularization term on the primal feasibility and thus may result in a lager primal residual of λ₀. Hence, an approach to adjust the penalty parameter r₀ and r₁ is to balance the residual of λ₀ is approximately equivalence to the residual of u.

2.6 Performance evaluations

Both quantitative and qualitative analysis will be included in our discussion for the evaluation of the TAC-TV algorithm. The root-mean-square error (RMSE), peak signal to noise ratio (PSNR), and normalized root mean square distance (NRMSD) are used as measures of the reconstruction quality. The RMSE, PSNR and NRMSD are defined as follows: $RMSE = \sqrt{\frac{\sum_{1}^{N} {| u_{i} - I_{i} |}^{2}}{N}},$ (32) $PNSR = 10 {log}_{10} (\frac{N {max}^{2} (I)}{\sum_{1}^{N} {| u_{i} - I_{i} |}^{2}})$ (33) $NRMSD = \sqrt{\frac{\sum_{1}^{N} {| u_{i} - I_{i} |}^{2}}{\sum_{1}^{N} I_{i}^{2}}}$ (34)

Where u and I denote the reconstruction image and the ideal phantom, respectively. And N denotes the total numbers of the pixels, and i is the index of the pixels. The smaller the RMSE and NRMSD value, the better the performance of image reconstruction. The higher PSNR shows that the reconstruction image is of higher quality.

In order to quantify the similar texture benefit of the proposed method, the often used Haralick texture distance (HTD) is taken [29]. The region of interest (ROI) is selected from the truth background images and the Euclidean distance between the textures of the truth background image and reconstructed images is used as the quantitative measure. The closer the distance is, the better the image texture preservation of the ROI region is.

3 Experimental studies

In this section, a series of experiments are conducted to investigate the proposed algorithm. It should be noted that as we have not proven and discussed the convergence of the TAC-TV algorithm, we first carry out the inverse-crime study to verify that for the pseudo real CT data the proposed method can reach the nearly converged and approximately exact reconstruction within the error tolerance. Moreover, the comparisons between the typical and the state-of-the-art methods are also conducted for both simulation dataset and the real CT projections.

The experiments are conducted on a PC equipped with dual-core Intel Core i5-2400 CPU at 3.10 GHz, with a memory size of 8.00 Gigabytes. The codes are written in MATLAB 2014b, and some of them are in C and CUDA with the “mex” interface provided by MATLAB and interacting with MATLAB code. The graphics processing unit (GPU) used is NVidia GeForce GTX 570 with 1.25 Gigabytes global memory and 480 CUDA cores.

In the running program, the choices of the values of some important parameters should be emphasized. In the proposed reconstruction model, the choices of a and b which control the relative strength of the first and the second order information of image. In our trials, it should be noted that it is usually a better choice to set a > b in order to achieve appropriate performances. In all of the experiments of this paper, we set a = 1 and b = 0.4.

3.1 Inverse-crime studies

In the subsection, the inverse-crime studies on the CT image of an anthropomorphic head phantom are performed to test the reconstruction precision and the convergence. We assume that if the reconstructed images are identical to the original CT image, the algorithm is seen as convergent. However, the exact reconstruction may take numerous iterations, and we argue that when the difference between the images is within the tolerance the reconstruction can be regarded as nearly convergent.

Here we present the original image in the left of Fig. 1. The projections are generated using a fast and parallel algorithm proposed by Gao [30], and thus the projection data b and the intended image u are noiseless and perfectly consistent. The scanning procedure is set as the standard fan beam geometry. 655 projections are acquired over 360 degrees along the standard circle trajectory. The distance of X-ray source to the iso-center is 1100 mm, and the distance of the center of the detector to the iso-center is 550 mm. The equivalent line detector has 1024 bins with size of 0.388 mm for each. Thus, the projection data has a data size of 1024×655. The reconstruction image has a size of 512×512 with each pixel having the size of 0.49 mm×0.49 mm. The generated projection data is as shown in the right of Fig. 1.

Fig.1

Background in truth image and the corresponding sinogram for inverse-crime studies.

In the inverse-crime studies, the stop criterion is set as: 1) if the change in RMSE is less than 10^-5 or 2) the iteration reach the maximum iteration number (set to 2000) the iteration is terminated. When either the above conditions is satisfied, the program terminated. The running time for each iterative loop is about 0.43 seconds. The reconstructed image is shown in Fig. 2, and the RMSE of the result compared to the standard image is 3.2E-4. For better visual comparison, both a full display window and a much narrower window are chosen for the result. The result suggests that for both windows the reconstructed image can be regarded as a nearly convergent result which implies that the proposed algorithm may reach a stratifying reconstruction.

Fig.2

Reconstruction and comparison of the inverse-crime studies: the top left is the reconstructed image in a display window of [0, 0.0629]. The top right is the error image between the reconstruction and the background truth image with the display window of [–0.0001, 0.0001]. The bottom left is the background in truth image in a much narrower window of [0.0095, 0.0285]. The bottom right is the reconstructed image in the same window as its left.

3.2 Comparisons of Sparse-view reconstructions

A series of tests under both noiseless and noisy datasets with different levels are conducted in this subsection. A slice of CT images from scanning of a rabbit head is taken into usage which is shown in Fig. 3. Four ROI regions are selected and marked with red rectangles in Fig. 3. We focus on the test of sparse-view reconstructions. The rotation radius is 164.86 mm, and the source to the center of the detector is 1052 mm. The phantom is of the size of 308×308 with 0.148 mm×0.148 mm for each pixel. The detector has 308 bins with length of 0.74 mm for each. The reconstruction programs are terminated when either of the following conditions is satisfied: 1) the residual for the feasibility conditions is less than 0.01 and 2) the iteration reaches the maximum number 2000. The reconstruction with noise free is shown in Fig. 4, which shows the residuals of the TAC-TV method is lower than 0.01 after about 1530 iterations.

Fig.3

Digital Rabbit head used in simulation study: from the left to right is truth background image in a display window of [0, 1] and [0, 0.34], the ROIs in display windows [0,0.34], respectively. The zooming images of ROIs are placed at the right portion.

Fig.4

Behavior of KKT residual norm for TAC-TV experiments from noise-free dataset.

3.2.1 Noise-free cases

In the first test, 72 views of projection, within the 360 degrees around the circle trajectory, are acquired for reconstruction. In the following, the performance of the proposed method by comparing it with the recent state-of-the-art algorithm TV-ADM is illustrated in Fig. 5. All methods could recover image well from sparse projection in visual inspection.

Fig.5

Image reconstruction (with zoomed ROIs) of the Rabbit head from noise-free projection dataset. Results of the A) TV-ADM, B) TAC-TV-ADM. Display window is [0, 0.34].

Table 3 lists the RMSE, PSNR, and the NRMSD measures of the image reconstructed by different algorithms with 2000 iterations. Since the reconstructed images in Fig. 5 have little difference from visual sense, the texture distance of ROI regions are employed and the results listed in Table 4 indicate that the proposed method outperforms TV method in image texture preservation of the ROI region.

Table 3

Evaluations of the results reconstructed by different method from noise-free projection data in digital Rabbit head experiments

	RMSE	PSNR	NRMSD
TV-ADM	9.459e-3	40.4834	6.73e-2
Proposed method	3.793e-3	48.4195	2.70e-2

Table 4

Texture distance of ROIs: Texture distance between the reference truth background image and the reconstructed images

	ROI 1	ROI 2	ROI 3	ROI 4
TV-ADM	0.7210	1.1876	1.3820	0.6653
Proposed method	0.3267	0.4536	0.5032	0.2696

To visualize the difference in detail, besides the texture distance of ROI regions, horizontal profiles of the resulting images (Fig. 6) are drawn along the 128th column from the 90th row to the 170th row. The images obtained by use of the proposed algorithm are reasonably accurate with only small distortions, and the zoom-in images may indicate that the proposed method can effectively capture the curvature of image.

Fig.6

Horizontal profiles (128th column) in the reconstruction results of the Rabbit head from noise-free projection dataset. Result of the A) TV-ADM and B) TAC-TV-ADM method.

3.2.2 Noisy cases

In this subsection, we carried out the experiments from the projection dataset at different noise levels to further illustrate the practicality of the algorithm. And the Poisson noise is considered: $prob (k, λ) = \frac{λ^{k}}{k!} e^{- λ}, λ = N_{0} exp (- p),$ (35) where k is the pixels in the projection data, N₀ denotes the incident X-ray intensity, and p is the normalized projection in real space. Here, the total view number of the experiment still was 72.

Three cases with different noise level are considered, and the initial numbers of photons N₀ are set to 3 × 10⁶, 10⁶, and 5 × 10⁵ for noise level 1, 2 and 3, respectively. For the noise level 1 and 2, parameter r₀ was set to 200, and for noise level 3, r₀ was set to 100 in TAC-TV-ADM algorithm. For the TAC-TV method, λ₁, λ₂, λ₃ and λ₄ are set to 1, 0.02, 180 and 180, respectively.

The images reconstructed by TV and TAC-TV methods from different noisy projection are shown in Fig. 7 and the RMSE, PSNR and NRMSD are listed in Table 5. The reconstructed images in Fig. 7 are very similar to each other, so the comparison of ROI regions are shown in Fig. 8 and the corresponding texture distance are listed in Table 6. Horizontal profiles of these images along the 128th column of the different noisy cases are shown in Figs. 9, 10 and 11, respectively. The presented TAC-TV image reconstruction has better performance on texture preservation, as expected.

Fig.7

Image reconstruction of the Rabbit head from noisy projection data. Images in the upper and the bottom row are the TV-ADM and the proposed method. From the left column to the right column are the reconstruction results from three different noise levels (level1, 2 and 3, respectively). Display window is [0, 0.34].

Fig.8

Reconstruction and comparison of the ROIs: from the top row to the bottom row are the TV-ADM reconstruction and the TAC-TV reconstruction by 2000 iterations. From the left column to the right column are the reconstruction results from three different noise levels (level 1, 2 and 3, respectively). Display window is [0, 0.34].

Fig.9

Horizontal profiles (128th column) in the reconstruction results of the Rabbit head from the noisy projection dataset (N₀ = 3 ×10⁶). Results of the A) TV-ADM and B) TAC-TV-ADM.

Fig.10

Horizontal profiles (128th column) in the reconstruction results of the Rabbit head from the noisy projection dataset (N₀ = 10⁶). Results of the A) TV-ADM and B) TAC-TV-ADM.

Fig.11

Horizontal profiles (128th column) in the reconstruction results of the Rabbit head from the noisy projection dataset (N₀ = 5 ×10⁵). Results of the A) TV-ADM and B) TAC-TV-ADM.

Table 5

Evaluations of the results reconstructed by different methods under three noise levels in digital Rabbit head studies

	Noise level 1			Noise level 2			Noise level 3
	RMSE	PSNR	NRMSD	RMSE	PSNR	NRMSD	RMSE	PSNR	NRMSD
TV-ADM	9.700e-3	40.2340	6.93e-2	9.900e-3	40.0896	7.04e-2	1.027e-2	39.7679	7.31e-2
TAC-TV-ADM	7.263e-3	42.7774	5.17e-2	8.633e-3	41.2771	6.14e-2	9.74e-3	40.3016	6.87e-2

Table 6

Texture distance of ROIs between the reference truth background image and the reconstructed images under three noise levels

		ROI 1	ROI 2	ROI 3	ROI 4
Noise 1	TV-ADM	0.7418	1.2205	1.4351	0.6803
	TAC-TV-ADM	0.5657	0.8709	1.0033	0.5223
Noise 2	TV-ADM	0.7458	1.2293	1.4447	0.6940
	TAC-TV-ADM	0.6729	1.0377	1.1668	0.6337
Noise 3	TV-ADM	0.7357	1.2826	1.4773	0.7571
	TAC-TV-ADM	0.7287	1.2196	1.3206	0.7295

3.3 Real CT data verifications

In this subsection, real CT verifications are carried out to test the proposed algorithm. These projections are acquired from an adult male rabbit in vivo on a CBCT imaging system. We declare and guarantee that this rabbit has never been tortured in the data acquisition and this data set is only for academic research. However, in order to obtain effective projection data, anesthetics with the minimum and safe dose is applied on this rabbit. The data acquisition takes only 4 minutes when this rabbit is under the quiet and rest state. The radiation dose and the exposure duration are strictly controlled to not make more harm to the rabbit. The breath of rabbit always existing, and thus we only acquired the projections of its head. More details of the data acquisition are omitted and we focus on the results of the experiments.

The distance of the source to the rotation center is 502.8 mm, and the distance of the source to the center of the detector is 1434.7 mm. The detector has 1024 bins with pixel size of 0.3810 mm for each. The reconstructed image is of 512×512 with 0.2671 mm×0.2671 mm for each image pixel. There are 360 projections acquired within the range of 360 degrees. Note that in this first test of the real CT data experiments all of the 360 projections are used for FBP, TV-ADM and the proposed method. For the second test, only 180 projections are taken to use with angle increment of 2 degrees.The reconstructions of FBP, TV-ADM and the proposed method from 360 views and 180 views are as shown in Figs. 12 and 13, respectively. We set the FDK result as the background reference image, and compared the RMSEs of the results of the iterative ones. In Fig. 9, the RMSEs for TV-ADM and the proposed method are 1.27E-4 and 0.94E-4. Moreover, in Fig. 10, the RMSEs are 2.7E-3 and 1.563E-3.

Fig.12

360 views reconstructions for FBP, TV-ADM and the proposed method: from the left column to the right column there are FBP reconstruction, TV-ADM reconstruction and the TAC-TV reconstruction by 500 iterations, respectively. The top and bottom row are the same image displayed windows [0, 0.133] and [0.0095, 0.0665].

Fig.13

180 views reconstructions for FBP, TV-ADM and the proposed method: from the left column to the right column there are FBP reconstruction, TV-ADM reconstruction and the TAC-TV reconstruction by 100 iterations, respectively. The top and bottom column are the same image displayed windows [0, 0.133] and [0.0095, 0.0665].

4 Discussion and conclusion

In this paper, we introduce a new TAC-TV model for CT image reconstruction for sparse-view data and develop an efficient and fast numerical algorithm based on the Augmented Lagrangian method and alternating direction method. The proposed algorithm was investigated by inverse-crime studies, sparse-view reconstructions experiments, and real CT data.

Different from the previous TV-based regularization strategies, the absolute value of curvature is introduced and incorporated in the TAC-TV model. The key motivation for the TAC-TV model is to utilize the second order information of the intended images for sparse-view reconstruction, which aims at reducing the radiation dose in sparse-view CT reconstruction and relieving the piecewise constant assumption of TV-based methods. Compared with the TV-based method, the proposed TAC-TV method shows better image quality in ROI regions and reducing the stair-case effect in sparse-view CT image reconstruction. In addition, the proposed method is robust to noise and shows better results than TV method. The view number of projection has significant influence on sparse-view CT reconstruction. Reconstruction under different views yields different levels of image quality. In this study, the view number is set as 72 and the gains from the proposed method are outstanding compared with those from TV method. By comparing the reconstructed image of ROI regions, the proposed method does have some advantages over the TV method. Parameter selection is another important issue in ADM scheme. Generally speaking, they are often empirically selected by quantity evaluation and visual inspection. In this paper, the parameter adapting strategy based on KKT conditions as a supplementary is investigated. The better selection strategy can help us achieve a good and stable performance of balancing the accuracy and efficiency.

Moreover, the convergence of the TAC-TV algorithm is verified by numerical experiments without strict theoretical proof. Future work will focus on the theoretical analysis of the proposed algorithm. Furthermore, more investigations for the comprehensive performances of the proposed algorithm should be studied as well.

Footnotes

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 61372172 and 61601518). We also want to thank the anonymous referees for giving us useful comments and suggestions to improve this paper. The contents of this paper are solely the responsibility of the authors and do not necessarily represent the official views of the NSFC.

References

Wernick

M.N.

and Aarsvold

J.N.

, Emission tomography: The fundamentals of PET and SPECT, 2004, (Academic Press).

Tuy

, An inversion formula for cone-beam reconstruction, SIAM J Apply Math 43 (1983), 546–552.

and Leahy

R.M.

, Iterative reconstruction techniques in emission computed tomography, Phys Med Biol 51 (2006), 541–578.

Donoho

, Compressed Sensing, IEEE Trans Inf Theory 52 (2006), 1289–1206.

Candes

, Romberg

and Tao

, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency Information, IEEE Trans Inf Theory 52 (2006), 489–509.

Candes

, Romberg

and Tao

, Stable signal recovery from incomplete and inaccurate measurements, Commun Pur Appl Math 59 (2006), 1207–1223.

Rudin

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992), 259–268.

Ertas

, Yildirim

, Kamasak

and Akan

, Iterative image reconstruction using non-local means with total variation from sufficient projection data, Journal of X-ray Science and Technology 24(1) (2016), 1–8.

Wang

, Yang

, Yin

and Zhang

, A new alternating minimization algorithm for total variation image reconstruction, SIAM Journal on Imaging Sciences 1(3) (2008), 248–272.

10.

Lysaker

, Lundervold

and Tai

X.C.

, Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time, Image Processing IEEE transactions 12(12) (2003), 1579–1590.

11.

Chan

R.H.

, Liang

, Wei

, Nikolova

and Tai

X.C.

, High-order total variation regularization approach for axially symmetric object tomography from a single radiograph, Inverse Problems and Imaging 9(1) (2015), 55–77.

12.

Chan

, Marquina

and Mulet

, High-order total variation-based image restoration, SIAM Journal on Scientific Computing 22(2) (2000), 503–516.

13.

Carlos

B.L.

and Chen

, Fast iterative algorithms for solving the minimization of curvature-related functionals in surface fairing, International Journal of Computer Mathematics 90(1) (2013), 92–108.

14.

Zhu

and Chan

, Image denoising using mean curvature of image surface, SIAM Journal on Imaging Sciences 5 (2012), 1–32.

15.

Carlos

B.L.

, Chen

and Victor

U.C.

, Image denoising using the Gaussian curvature of the image surface, Numerical Methods for Partial Differential Equations 32(3) (2015).

16.

Chumchob

, Chen

and Loeza

C.B.

, A fourth-order variational image registration model and its fast multi-grid algorithm, SIAM J Multiscale Model Simul 9(1) (2011), 89–128.

17.

Tasdizen

, Whitaker

, Burchard

and Osher

, Geometric surface processing via normal maps, ACM Trans Graph 22(4) (2003), 1012–1033.

18.

Loeza

C.B.

and Chen

, Fast numerical algorithms for Euler-elastica inpainting model, Int J Modern Math 5(2) (2010), 157–182.

19.

Yang

, Chen

and Yu

, Homotopy method for a mean curvature-based denoising model, Appl Numer Math 62(3) (2012), 185–200.

20.

Tai

X.C.

, Hahn

and Chung

G.J.

, A fast algorithm for Euler’s Elastica model using augmented Lagrangian method, SIAM J Imaging Sciences 4(1) (2011), 313–344.

21.

Zhu

, Chan

and Tai

X.C.

, Augmented Lagrangian method for a mean curvature based image denoising model, Inverse Problems and Imaging 7 (2013), 1409–1432.

22.

Zhang

, Wang

and Yan

, Image reconstruction based on total-variation minimization and alternating direction method in linear scan computed tomography, Chin Phys B 22(7) (2013), 582–589.

23.

Wang

, Zhang

and Cai

, Image reconstruction algorithm based on inexact alternating direction total-variation minimization, Acta Phys Sin 62 (2013) 198701, 1–6.

24.

Chan

T.F.

and Shen

, Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods, SIAM, Philadelphia, 2005.

25.

Mumford

, Elastica and computer vision. In Algebraic Geometry and Its Applications, Bajaj

C.L.

, ed., Springer-Verlag, New York (1994), 491–506.

26.

, Yin

and Jiang

, An efficient augmented Lagrangian method with applications to total variation minimization, Computational Optimization and Applications 56 (2013), 507–530.

27.

Lions

P.L.

and Mercier

, Splitting algorithms for the sum of two nonlinear operators, SIAM J Numer Anal 16 (1979), 964–979.

28.

Boyd

, Parikh

, Chu

, Peleato

and Eckstein

, Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine Learning 3(1) (2011), 1–122. [Online]. Available: https://dx-doi-org.web.bisu.edu.cn/10.1561/2200000016

29.

Zhang

, Han

, Liang

Z.R.

, et al. Extracting information from previous full-dose CT scan for knowledge-based Bayesian reconstruction of current low-dose CT images, IEEE Trans Medical Imaging 35(3) (2016), 860–870.

30.

Gao

, Fast parallel algorithm for the x-ray transform and its adjoint, Med Phys 39 (2012), 7110–7120.