Ordered-subset Split-Bregman algorithm for interior tomography

Abstract

Inspired by the Compressed Sensing (CS) theory, it has been proved that the interior problem of computed tomography (CT) can be accurately and stably solved if a region-of-interest (ROI) is piecewise constant or polynomial, resulting in the CS-based interior tomography. The key is to minimize the total variation (TV) of the ROI under the constraint of the truncated projections. Coincidentally, the Split-Bregman (SB) method has attracted a major attention to solve the TV minimization problem for CT image reconstruction. In this paper, we apply the SB approach to reconstruct an ROI for the CS-based interior tomography assuming a piecewise constant imaging model. Furthermore, the ordered subsets (OS) technique is used to accelerate the convergence of SB algorithm, leading to a new OS-SB algorithm for interior tomography. The conventional OS simultaneous algebraic reconstruction technique (OS-SART) and soft-threshold filtering (STF) based OS-SART are also implemented as references to evaluate the performance of the proposed OS-SB algorithm for interior tomography. Both numerical simulations and clinical applications are performed and the results confirm the advantages of the proposed OS-SB method.

Keywords

Ordered subset Split-Bregman interior tomography compressive sensing total variation minimization piecewise constant imaging model

1 Introduction

In the computed tomography (CT) field, the interior problem refers to reconstruct a region of interest (ROI) only from data associated with x-rays through the ROI. It has a great potential to handle large objects, minimize radiation dose, reduce system cost, enhance temporal resolution and increase scanner throughput, etc. Because the interior problem does not have a unique solution in an unconstrained setting, an internal ROI cannot be accurately reconstructed only from the truncated data based on the conventional CT theory. In 2007, it was proved and demonstrated that the interior problem can be exactly and stably solved if a sub-region inside the ROI is known [1 –4], which is referred to as knowledge-based interior tomography. However, it is very difficult to obtain precise knowledge of a sub-region in many practical applications, such as cardiac imaging. Inspired by the elegant compressive sensing (CS) theory [5, 6], in 2009 Yu et al. proved that exact interior reconstruction is achievable by minimizing the total variation (TV) of an ROI if the ROI is piecewise constant [7] or polynomial [8], leading to the CS-based interior tomography. Those theoretical results were also verified and evaluated using clinical projections [9].

The main idea of CS is that most signals or images are sparse in an appropriate representation, that is, a majority of their coefficients are close or equal to zero. It is well known that the X-ray attenuation coefficient often varies mildly within an anatomical component and jumps significantly around the borders of tissue structures. The discrete gradient transform (DGT), whose L1-norm is referred to as TV, has been widely utilized as a sparsifying operator, and the corresponding TV minimization method has been adopted for CS-based CT image reconstruction including interior tomography. The TV minimization can be implemented by the steepest descent (SD) method or other optimization methods[10 –12]. In 2010, Yu and Wang applied a soft-threshold filtering (STF) framework to minimize TV by constructing a pseudo-inverse transform for DGT [13].

Coincidentally, the Split-Bregman (SB) approach has received more attentions to solve a broad class of L1-regularized optimization problems [14]. This method applies the Bregman update [12] after introducing an auxiliary variable to the original L1-regularized problem. As a result, it ‘splits’ the L1 and L2 components of the cost function into several sub-problems. This method converges very fast if the sub-problems can be efficiently solved with proper parameters. However, when it is applied to the CT reconstruction problem, several forward/back-projection operations are required in each iteration, which are very time-consuming [15 –18]. It is well known that the ordered-subset (OS) method is effective to accelerate the convergence speed of the image reconstruction by grouping the projections into some ordered subsets [19 –21]. Motivated by the algorithmic application of the CS-based interior tomography and the merits of SB and OS techniques, in this paper we develop an ordered-subsets SB (OS-SB) method and apply it for CS-based interior reconstruction and few-view reconstruction. During the review process, one of the reviewers brought our attention to the work by Zhao et al. [22]. Although the goal of our work is similar to that of [22], there are several significant differences. First, the object function of [22] is constructed based on a modified high-order TV while ours is based on the ordinary TV. Second, the reference [22] applied the OS simultaneous algebraic reconstruction technique (OS-SART) and SB techniques to solve the unknown variables in order in the outer iteration, while our method applies the OS technique directly to SB and it results in the OS-SB. Third, the OS-SART and STF based OS-SART are also implemented as references to evaluate the performance of the proposed OS-SB algorithm for interior tomography.

The rest of this paper is organized as follows. In Section 2, we present the CT imaging model, TV minimization problem and review the STF algorithm. The SB and OS-SB methods for CT image reconstruction with TV minimization are introduced in Section 3. In Section 4, we perform numerical simulations and clinical applications to evaluate the proposed OS-SB algorithm in reference to the conventional OS-SART and STF based OS-SART. Finally, we conclude this paper in the last section.

2 Imaging model and related work

2.1 Image reconstruction model

The image reconstruction problem in CT can be modeled by the following linear equation: $Af = y + e,$ (1) where $f = (f_{j})_{J \times 1} \in ℝ^{J}$ is an image vector to be reconstructed, $y = (y_{i})_{I \times 1} \in ℝ^{I}$ is a measured projection data vector, $A = (a_{ij})_{I \times J} \in ℝ^{I \times J}$ is a projection system matrix and $e = (e_{i})_{I \times 1} \in ℝ^{I}$ is a measurement noise vector. For an accurate fan-beam geometry using an area integral model [23], the j-thpixel can be viewed as a rectangular region with a constant value f_j, the i-th measured datum y_i is an integral of the weighted product between the x-ray linear attenuation coefficient and the corresponding area of pixel partially covered by the narrow fan-beam from the x-ray source to the detector element. Therefore, the component a_ij in Equation (1) can be viewed as a normalized interaction area between the j-th pixel and i-th fan-beam path. The goal of CT image reconstruction is to recover the image f from the noisy data y + e .

The representative SART algorithm [24] can be used to solve (1) and it can be written as ${\hat{f}}_{j}^{(k + 1)} = {\hat{f}}_{j}^{(k)} + λ_{k} \cdot \frac{1}{a_{+ j}} \sum_{i = 1}^{I} \frac{y_{i} - \sum_{j = 1}^{J} a_{ij} {\hat{f}}_{j}^{(k)}}{a_{i +}} \cdot a_{ij},$ (2) where k is the iteration index, i = k mod (I) +1, $a_{+ j} = \sum_{i = 1}^{I} a_{ij}$ , $a_{i +} = \sum_{j = 1}^{J} a_{ij}$ , and 0 < λ_k < 2 is a free relaxation parameter. Let $Λ^{+ J} \in ℝ^{J} \times ℝ^{J}$ be a diagonal matrix with $Λ_{jj}^{+ J} = \frac{1}{a_{+ j}}$ and $Λ^{I +} \in ℝ^{I} \times ℝ^{I}$ be a diagonal matrix with $Λ_{ii}^{I +} = \frac{1}{a_{i +}}$ , then (2) can be simplified as ${\hat{f}}^{k + 1} = {\hat{f}}^{k} + λ_{k} (Λ^{+ J} A^{T} Λ^{I +} (y - A {\hat{f}}^{k})) .$ (3)

In (3), the weighting matrices Λ ^+J and Λ ^I+ play a role to normalize the matrix A ^T A . In fact, we can ignore the weighting matrices Λ ^+J and Λ ^I+ and re-express (3) as ${\hat{f}}^{k + 1} = {\hat{f}}^{k} + λ_{k} A^{T} (y - A {\hat{f}}^{k}) .$ (4)

Eq.(4) is the well-known steepest gradient descent algorithm, and λ_k plays a role to control the length of numerical search step whose value should be empirically selected.

2.2 DGT and TV minimization

Although most of the medical images are not naturally sparse in spatial domain, they can be sparsely expressed in an appropriate domain by applying certain transform. One of the most notable sparse transforms is DGT. It is well known that a two dimensional image vector f = (f_m,n) _M×N can be re-arranged into a column vector $f = (f_{j})_{J \times 1} \in ℝ^{J}$ with the relationship $f_{j} = f_{m, n}, j = (m - 1) \times N + n, 1 \leq m \leq M, 1 \leq n \leq N, J = M \times N .$ (5)

In this paper, both (f_m,n) _M×N and (f_j) _J×1 are used to represent an image for convenience. Here, we also assume that an image (f_m,n) _M×N satisfies the so-called Neumann condition on the boundary: $\begin{matrix} f_{0, n} = f_{1, n} and f_{M, n} = f_{M + 1, n} for 1 \leq n \leq N, \end{matrix}$ $\begin{matrix} f_{m, 0} = f_{m, 1} and f_{m, N} = f_{m, N + 1} for 1 \leq m \leq M . \end{matrix}$

Then the DGT of a two-dimensional image f = (f_m,n) _M×N can be defined as

$\begin{matrix} \nabla f_{m, n} & = & \sqrt{(\nabla_{x} f)_{m, n}^{2} + (\nabla_{y} f)_{m, n}^{2}} \\ = & \sqrt{{(\frac{f_{m, n} - f_{m + 1, n}}{δ_{x}})}^{2} + {(\frac{f_{m, n} - f_{m, n + 1}}{δ_{y}})}^{2}}, 1 \leq m \leq M, 1 \leq n \leq N, \end{matrix}$ (6) and the corresponding TV can be defined as ${∥ f ∥}_{TV} = | \nabla f | = \sum_{m = 1}^{M} \sum_{n = 1}^{N} | \nabla f_{m, n} |,$ (7) where δ_x and δ_y represent the sampling intervals. In practice applications, usually δ_x = δ_y and they can be absorbed by other parameters. The TV minimization based CT image reconstruction can be expressed as a constrained optimization problem $min {∥ f ∥}_{TV} s . t ∥ Af - y ∥ \leq ɛ,$ (8) where ɛ is an estimated noise level.

2.3 Soft-threshold filtering algorithm

The problem (8) is often solved by the alternating minimization scheme which includes two major steps. The first step is to impose the projection data constrain by using an iteration method such as algebraic reconstruction technique (ART) [25], SART, expectation maximization (EM) [26], etc. The second step is to minimize the image’s TV by the SD method, STF method, etc. Therefore, there are many algorithms to solve (8) by combining different components in the two steps. In this paper, the STF based SART algorithm is chosen as the reference. For completeness, the STF technique for TV regularization is summarized as follows.

The convergence of the STF approach [27] was originally proved and applied for the L1 norm regularized minimization of an invertible sparse transform. However, the DGT is not invertible and it does not satisfy the restricted isometry property (RIP) required by the CS theory. In other words, the STF algorithm cannot be directly applied for TV regularized CT image reconstruction. Motivated by this challenging situation, Yu and Wang constructed several pseudo-inverse transforms and applied the STF approach to image reconstruction from a limited number of projections [13]. This method was also adopted for interior tomography [9]. Assume that ${\hat{f}}^{k}$ is the updated image after the iteration of the SART, then we can obtain $\nabla {\hat{f}}_{m, n}^{k} = \sqrt{{({\hat{f}}_{m, n}^{k} - {\hat{f}}_{m + 1, n}^{k})}^{2} + {({\hat{f}}_{m, n}^{k} - {\hat{f}}_{m, n + 1}^{k})}^{2}} .$ (9)

One possible pseudo-inverse transform for the DGT can be constructed as follow ${\tilde{f}}_{m, n}^{k} = \frac{1}{4} (2 {\tilde{f}}_{m, n}^{k, a} + {\tilde{f}}_{m, n + 1}^{k, a} + {\tilde{f}}_{m + 1, n}^{k, a}),$ (10) ${\tilde{f}}_{m, n}^{k, a} = {\begin{matrix} \frac{2 {\hat{f}}_{m, n}^{k} + {\hat{f}}_{m + 1, n}^{k} + {\hat{f}}_{m, n + 1}^{k}}{4} \nabla {\hat{f}}_{m, n}^{k} < ω \\ {\hat{f}}_{m, n}^{k} - \frac{ω (2 {\hat{f}}_{m, n}^{k} - {\hat{f}}_{m + 1, n}^{k} - {\hat{f}}_{m, n + 1}^{k})}{4 \nabla {\hat{f}}_{m, n}^{k}} \nabla {\hat{f}}_{m, n}^{k} \geq ω \end{matrix},$ (11) ${\tilde{f}}_{m + 1, n}^{k, a} = {\begin{matrix} \frac{{\hat{f}}_{m + 1, n}^{k} + {\hat{f}}_{m, n}^{k}}{2} \nabla {\hat{f}}_{m, n}^{k} < ω \\ {\hat{f}}_{m + 1, n}^{k} - \frac{ω ({\hat{f}}_{m + 1, n}^{k} - {\hat{f}}_{m, n}^{k})}{2 \nabla {\hat{f}}_{m, n}^{k}} \nabla {\hat{f}}_{m, n}^{k} \geq ω \end{matrix},$ (12) ${\tilde{f}}_{m, n + 1}^{k, a} = {\begin{matrix} \frac{{\hat{f}}_{m, n + 1}^{k} + {\hat{f}}_{m, n}^{k}}{2} \nabla {\hat{f}}_{m, n}^{k} < ω \\ {\tilde{f}}_{m, n + 1}^{k} - \frac{ω ({\hat{f}}_{m, n + 1}^{k} - {\hat{f}}_{m, n}^{k})}{2 \nabla {\hat{f}}_{m, n}^{k}} \nabla {\hat{f}}_{m, n}^{k} < ω \end{matrix},$ (13) where ω is a threshold for STF, and the sampling interval parameter δ_x = δ_y is absorbed by ω.

3 Ordered subsets Split-Bregman

3.1 Split-Bregman method

To deal with the difficulty of solving the constrained optimization problem, (8) is usually converted into an approximately equivalent unconstrained convex minimization problem using a quadratic penalty function $min {∥ f ∥}_{TV} + \frac{μ}{2} {∥ Af - y ∥}_{2}^{2},$ (14) where μ is a constant parameter to balance the two terms in (14). The SB method uses a variable-splitting technique as well as the Bregman iteration to solve the non-differentiable regularizers such as the TV minimization problem. The major steps of the SB for CT reconstruction problem (8) can be summarized as follows. Two new variables are firstly introduced to convert (14) to an equivalent constrained convex minimization problem $min_{f, d_{x}, d_{y}} {∥ (d_{x}, d_{y}) ∥}_{2} + \frac{μ}{2} {∥ Af - y ∥}_{2}^{2} s . t d_{x} = \nabla_{x} f, d_{y} = \nabla_{y} f .$ (15) where ${∥ (d_{x}, d_{y}) ∥}_{2} = \sum_{m, n} \sqrt{d_{x, m, n}^{2} + d_{y, m, n}^{2}}$ .

Problem (15) can be then converted into another unconstrained optimization by adding penalty function terms into (14) $min_{f, d_{x}, d_{y}} {∥ (d_{x}, d_{y}) ∥}_{2} + \frac{μ}{2} {∥ Af - y ∥}_{2}^{2} + \frac{λ}{2} {∥ d_{x} - \nabla_{x} f ∥}_{2}^{2} + \frac{λ}{2} {∥ d_{y} - \nabla_{y} f ∥}_{2}^{2},$ (16) where λ plays the same role as μ. The problem (16) can be solved by using the Bregman iteration[12]

$\begin{matrix} (f^{k + 1}, d_{x}^{k + 1}, d_{y}^{k + 1}) = min_{f, d_{x}, d_{y}} {∥ (d_{x}, d_{y}) ∥}_{2} + \frac{μ}{2} {∥ A f - y ∥}_{2}^{2} + \frac{λ}{2} {∥ d_{x} - \nabla_{x} f - b_{x}^{k} ∥}_{2}^{2} \\ + \frac{λ}{2} {∥ d_{y} - \nabla_{y} f - b_{y}^{k} ∥}_{2}^{2}, \end{matrix}$ (17) $b_{x}^{k + 1} = b_{x}^{k} + (\nabla_{x} f^{k + 1} - d_{x}^{k + 1}),$ (18) $b_{y}^{k + 1} = b_{y}^{k} + (\nabla_{y} f^{k + 1} - d_{y}^{k + 1}),$ (19) where k represents the iteration index. For the above minimization problem (17), it can be solved by alternately performing iteration to update f and d _x and d _y $\begin{matrix} Step 1 f^{k + 1} = \underset{f}{arg min} \frac{μ}{2} {∥ Af - y ∥}_{2}^{2} + \frac{λ}{2} {∥ d_{x}^{k} - \nabla_{x} f - b_{x}^{k} ∥}_{2}^{2} + \frac{λ}{2} {∥ d_{y}^{k} - \nabla_{y} f - b_{y}^{k} ∥}_{2}^{2}, \end{matrix}$ (20) $\begin{matrix} Step 2 (d_{x}^{k + 1}, d_{y}^{k + 1}) = arg min_{d_{x}, d_{y}} {∥ (d_{x}, d_{y}) ∥}_{2} + \frac{λ}{2} {∥ d_{x} - \nabla_{x} f^{k + 1} - b_{x}^{k} ∥}_{2}^{2} + \frac{λ}{2} {∥ d_{y} - \nabla_{y} f^{k + 1} - b_{y}^{k} ∥}_{2}^{2} . \end{matrix}$ (21)

Obviously, the effectiveness of the alternating minimization approach depends on how fast the two subproblems (20) and (21) be solved.

In the step1 of (20), because the cost function is differentiable, various optimization techniques can be used, such as SD, Gauss-Seidel and conjugate gradient (CG) methods. In the step2 of (21), an analytic solution is available by using a generalized shrinkage formula [28] $d_{x}^{k + 1} = max (s^{k} - \frac{1}{λ}, 0) \frac{\nabla_{x} f^{k} + b_{x}^{k}}{s^{k}},$ (22) $d_{y}^{k + 1} = max (s^{k} - \frac{1}{λ}, 0) \frac{\nabla_{y} f^{k} + b_{y}^{k}}{s^{k}},$ (23) where $s^{k} = \sqrt{{| \nabla_{x} f^{k} + b_{x}^{k} |}^{2} + {| \nabla_{y} f^{k} + b_{y}^{k} |}^{2}} .$ (24)

The choice of the parameters μ and λ plays an important role in the implementation of the algorithm [14 , 29]. In (14), the parameter μ balances the data fidelity and TV term. If μ is too large, the reconstructed image is noisy. On the other hand, if μ is too small, the reconstructed image is blocky, and some finer structures may be lost which may provide critical information for clinical diagnosis. The parameter λ determines the convergence to the solution. If λ is too small, the corresponding penalty functions will be violated excessively. On the other hand, if λ is too large, the amount of update applied to d _x and d _y (and subsequently to f ) is small, and it implies tiny iterative steps and more iterations for convergence. Therefore, the parameter λ should be chosen in a way that balances the convergence speed and reconstruction accuracy. Unfortunately, there is no analytic method to optimize the parameters for iterative reconstruction. Following the classical paradigms in iterative CT reconstruction, in our experiments the parameters μ and λ are empirically selected based on extensive tests.

3.2 OS-SB method

Although the SB method converges quickly with effective solvers for the sub-problems (20) and (21), it is still slow due to the iterative updates which require several time-consuming forward/back-projection operations. The OS method is very useful to accelerate a reconstruction algorithm which employs all the measurements and updates the results simultaneously at each iteration [20]. The basic idea of OS-type algorithms is to group the measurements into an ordered sequence of subsets and utilize only one subset for each update/subiteration. Because the SB algorithm uses all the measurements to solve the subproblem (20), it is relatively straight forward to modify the SB algorithm to incorporate the OS idea. The resulting algorithm can be called as OS-SB method.

The proposed OS-SB method first groups the measured projections index {1, 2, … , I } into T ordered subsets {S₁, … , S_T }. Then, the standard SB algorithm is applied to each of the subsets in order which uses the rows of the projection matrix corresponding to these ray sums. An update using a single subset is called a subiteration. The resulting reconstruction result becomes the starting value for the next subiteration with respect to the next subset. One iteration of the OS-SB is defined as a single pass-through for all the specified subsets. Multiple iterations can be performed by passing through the same ordered subsets using as a starting point provided by the previous iteration. Considering the particularity of matrix A for CT reconstruction, in this paper we adopt the SD method [30] to solve (20). The detailed algorithm is summarized as the following pseudo-codes:

Initialization: ${\hat{f}}^{0} = 0, d_{x}^{0} = d_{y}^{0} = b_{x}^{0} = b_{y}^{0} = 0, k = 1$ ;

While the convergent criteria are not satisfied ${\hat{f}}^{k, 0} \leftarrow {\hat{f}}^{k - 1}, d_{x}^{k, 0} \leftarrow d_{x}^{k - 1}, d_{y}^{k, 0} \leftarrow d_{y}^{k - 1}, b_{x}^{k, 0} \leftarrow b_{x}^{k - 1};$

For t = 1, 2, …, T (For each subset S₁, …, S_T , updating the image by SB) ${\hat{f}}^{k, t} \leftarrow {\hat{f}}^{k, t - 1}, d_{x}^{k, t} \leftarrow d_{x}^{k, t - 1}, d_{y}^{k, t} \leftarrow d_{y}^{k, t - 1}, b_{x}^{k, t} \leftarrow b_{x}^{k, t - 1};$

Computing $s^{k, t} = \sqrt{| \nabla_{x} {\hat{f}}^{k, t} + b_{x}^{k, t} | + | \nabla_{y} {\hat{f}}^{k, t} + b_{y}^{k, t} |}$ ; ${\hat{f}}^{k, t, 0} \leftarrow {\hat{f}}^{k, t};$ For n = 1 to N (Implementing the SD method) ${\hat{f}}^{k, t, n} \leftarrow {\hat{f}}^{k, t, n - 1};$ $\begin{matrix} g^{n} & = & μ \cdot A_{S_{t}}^{T} (A_{S_{t}} {\hat{f}}^{k, t, n} - y_{S_{t}}) - λ \cdot \nabla_{x}^{T} (d_{x}^{k, t} - \nabla_{x} {\hat{f}}^{k, t, n} - b_{x}^{k, t}) \\ - λ \cdot \nabla_{y}^{T} (d_{x}^{k, t} - \nabla_{y} {\hat{f}}^{k, t, n} - b_{y}^{k, t}); \\ α_{n} & = & \frac{g^{n'} \cdot g^{n}}{g^{n'} \cdot \nabla_{f} g^{n} \cdot g^{n}}; \end{matrix}$ ${\hat{f}}^{k, t, n} = {\hat{f}}^{k, t, n} - α_{n} \cdot g^{n};$

End(n); $\begin{matrix} {\hat{f}}^{k, t} \leftarrow {\hat{f}}^{k, t, N}; \\ d_{x}^{k, t} = max (s^{k, t} - \frac{1}{λ}, 0) \frac{\nabla_{x} {\hat{f}}^{k, t} + b_{x}^{k, t}}{s^{k, t}}; \\ d_{y}^{k, t} = max (s^{k, t} - \frac{1}{λ}, 0) \frac{\nabla_{y} {\hat{f}}^{k, t} + b_{y}^{k, t}}{s^{k, t}}; \\ b_{x}^{k, t} = b_{x}^{k, t} + (\nabla_{x} {\hat{f}}^{k, t} - d_{x}^{k, t}); \\ b_{y}^{k, t} = b_{y}^{k, t} + (\nabla_{y} {\hat{f}}^{k, t} - d_{y}^{k, t}); \end{matrix}$ End(subsets); ${\hat{f}}^{k} \leftarrow {\hat{f}}^{k, T};$ $k \leftarrow k + 1;$

End (while).

In the above pseudo-code, y _{S
_t} represents the projection vector for the t-th subset and A _{S
_t} is corresponding system matrix. The proposed OS-SB method degenerates to the SB algorithm when T = 1. Because the projections and backprojections are performed only for the elements within a single subset, the computation cost for each subset in an OS-SB algorithm with T subsets roughly takes 1/T time of that treats the whole dataset as one subset. It also implies that the image can be updated T times per iteration.

3.3 Competitive method

To evaluate the performance of the SB and OS-SB for interior tomography, the SART and STF based SART are chosen for comparison in the experiments. Because the convergence of SART is slow, we also incorporate the OS technique into the SART approach resulting in an OS-SART [31]. It can be written as ${\hat{f}}_{j}^{(k + 1)} = {\hat{f}}_{j}^{(k)} + λ_{k} \cdot \sum_{i \in S_{t}} \frac{a_{ij}}{\sum_{i^{'} \in S_{t}} a_{i^{'} j}} \frac{y_{i} - \sum_{j = 1}^{J} a_{ij} {\hat{f}}_{j}^{(k)}}{a_{i +}},$ (25) where S_t represents the index set for the t-th subset, ⋃_1≤t≤TS_t = {1, 2, …, I}, t = k mod (T) +1. Equation (25) also has the corresponding matrix form ${\hat{f}}^{k + 1} = {\hat{f}}^{k} + λ_{k} (Λ_{S_{t}}^{+ J} A_{S_{t}}^{T} Λ_{S_{t}}^{I +} (y_{S_{t}} - A_{S_{t}} {\hat{f}}^{k})) .$ (26)

The STF based OS-SART (STF-OS-SART) can be summarized in the following pseudo codes:

Initialization: ${\hat{f}}^{0} = 0$ , k = 0 ; While the stopping criteria are not met ${\hat{f}}^{k, 0} \leftarrow {\hat{f}}^{k - 1}$

For t = 1, 2, …, T (For each subset S₁, …, S_T updating the image by the OS-SART) ${\hat{f}}^{k, t} = {\hat{f}}^{k, t - 1} + λ_{k} (Λ_{S_{t}}^{+ J} A_{S_{t}}^{T} Λ_{S_{t}}^{I +} (y_{S_{t}} - A_{S_{t}} {\hat{f}}^{k, t - 1}));$

End (subsets); ${\hat{f}}^{k} \leftarrow {\hat{f}}^{k, T};$

Performing the soft-threshold filtering for ${\hat{f}}^{k}$ using Equation (10) $k \leftarrow k + 1;$ End (while).

For the algorithms OS-SART, STF-OS-SART and OS-SB, the OS-SART is to only solve (1) without any regularization, and the STF-OS-SART and OS-SB are to solve the TV minimization problem (8). The STF-OS-SART employs an alternating minimization scheme. Because the STF does not work until the gradient is greater than the threshold, there are no differences between the STF-OS-SART and OS-SART in early iterations. The OS-SB converts the constrained problem (8) into an unconstrained problem (16) and solves it using the Bregman iteration. The major difference between the STF-OS-SART and OS-SB is that the STF-OS-SART performs STF after the image update using the OS-SART, while the OS-SB updates the image under the penalty function for constraints.

4 Experiment and result

4.1 Experiment design

The major goal of this paper is to evaluate the performance of the SB and OS-SB for interior tomography. Meanwhile, the sparse-view CT reconstruction is also considered. The OS-SART and the STF-OS-SART are implemented for comparison. For brevity, we denote the OS-SB algorithm using n subsets by OS-SB-n, the OS-SART algorithm using n subsets by OS-SART-n, and the STF based OS-SART-n as STF-OS-SART-n. All the algorithms are implemented in a hybrid mode of Matlab and C++. While the interface is implemented in Matlab, all the extensive computational parts are implemented in C++and complied via MEX function.

The mean square error (MSE) and structural similarity (SSIM) [32] indexes are used to quantitatively evaluate the performance of the proposed algorithm. The MSE is used to estimate perceived errors by averaging the squared intensity differences of reconstructed and reference image pixels. The SSIM quantitatively characterizes image degradation as perceived change in structural information and can be considered as a complementary measure to the MSE. The closer to 1.0 the value of SSIM is, the more fine structures are reconstructed.

4.2 Numerical simulation

A modified Shepp-Logan phantom whose parameters are from the source codes of Matlab is used for extensive numerical simulations. The phantom includes a series of piecewise constant ellipses and consists of 256×256 pixels, each of which covers an area of 0.781×0.781 mm². A circular scanning trajectory of radius 500 mm is assumed for typical fan-beam geometry. An equi-angular detector including 320 elements is also assumed with a fan-angle of 0.49 radian. For a full scan, different numbers of noise-free and noisy projections are equi-angularly collected to evaluate the proposed algorithm. To evaluate the OS-SB based interior reconstruction technique, each projection is truncated by discarding 80 detector elements on two sides of the sinogram, and the truncated projections covers a disk-shaped ROI with a diameter of 121.83 mm as shown in Fig. 1.

To test the OS-SB performance for CS-based interior tomography, 550 noise-free projections are uniformly acquired assuming a full scan. The images are reconstructed by the OS-SB assuming 1, 5 and 10 subsets, respectively. For the OS-SART and STF-OS-SART, the case of 10 subsets is selected for a fairly comparison. The parameters choosing for different algorithms are as following. The parameter λ_k in the OS-SART and STF-OS-SART is chosen as a constant 1.0. In the STF-OS-SART, an optimal threshold ω is determined by the projected gradient method for each filtering [13]. For the OS-SB, the parameter μ and λ are related to many factors such as the subset number, the projection scale, noise level, etc. In this study, the parameter λ is chosen as 10 for all numerical simulation experiments. With the increase of the subset number of the OS-SB, the projection number in each subset decreases. In order to balance the data fidelity and the regularization term, we have to choose a greater μ for a larger subset number. Regarding the discrete gradient, it is computed using Equation (9) and the sampling interval parameter δ_x = δ_y is absorbed by λ and μ. The final values of the parameters λ and μ of the OS-SB are summarized as in Table 1. The reconstructed results after 30 iterations are shown in Fig. 2 in a display window of [0.1 0.5]. Figure 3 plots the MSE and SSIM convergent curves with respect to the iteration number for different algorithms. Comparing with the OS-SART, the STF-OS-SART guarantees an accurate reconstruction as expected although there is almost no difference in earlier iterations. The OS-SB algorithm shows superiority at the beginning. The greater the subset number is used, the faster the convergent rate is achieved. It also demonstrates that the SB method with proper parameters is very fast and the OS-SB can further accelerate the convergence rate. In order to further compare the accuracy of the reconstruction, the profiles along the central horizontal and vertical lines of the reconstructed images are shown in Fig. 4 (here we only list the results of STF-OS-SART and OS-SB algorithms). Figure 4 clearly shows that the OS-SB algorithms outperform the STF-OS-SART-10 in reference to the truth of the original Shepp-Logan phantom. It is also observed that increasing the subset number can introduce noise in the reconstructed images because the projection number in each subset of the OS-SB with a greater subset number (OS-SB-5, OS-SB-10) is significantly smaller than that with a smaller subset number(OS-SB-1).

To test the feasibility of the OS-SB approach for sparse-view interior reconstruction, 60 views are uniformly sampled from a full scan range. Images are reconstructed by the OS-SB assuming 1, 2 and 5 subsets, respectively. The case of 5 subsets is selected for the OS-SART and STF-OS-SART to compare. The reconstructed results after 50 iterations are shown in Figs. 5–7 which can draw the similar conclusions as those reconstructed from 550 projections. However, the SSIM of OS-SART-5 diverges gradually with the increase of the iteration number, and the SSIM of OS-SB-5 decreases at the 49-th iteration. The reasons are as follow. First, the SART without the TV constrain cannot guarantee an exact and stable solution for interior problem especially for sparse-view case. Second, the OS algorithm does not converge to the true optimum for a subset number greater than one. Third, the MSE measures the L2 norm between the intermediate results and the ground truth f ^* as ${∥ {\hat{f}}^{(k)} - f^{*} ∥}_{2}^{2}$ , and the objective function for image reconstruction in (14) includes the TV term and data fiedlity term, while

${∥ A {\hat{f}}^{(k)} - y ∥}_{2}^{2} = {∥ A {\hat{f}}^{(k)} - A f^{*} ∥}_{2}^{2} \leq {∥ A ∥}_{2}^{2} {∥ {\hat{f}}^{(k)} - f^{*} ∥}_{2}^{2} .$ (27)

Therefore, we can see that MSE curve should be close but not necessary consistant with the convergence curve of the algorithm. Meanwhile, the SSIM describes the structure similarity of the images, and it not necessary consistent with the convergence curve of the algorihtm, either.

In practical applications, measurement noise is unavoidable. To test the feasibility of the proposed algorithm against data noise, the aforementioned tests are repeated for the Shepp-Logan phantom with 60 projection views corrupted by Gauss noise whose mean is 0 and standard deviation is 0.005. The PSNR for the projections is 42.08 dB. This noise model can be viewed as the special case of a dedicated bowtie filtering that the expected photon number arriving at each detector cells are the same. That is, each detector cell expects to receive the same number of photons. In order to decrease the influence of noise, in the experiment we choose a smaller value for μ and remain λ = 10 compared with the reconstruction from noiseless projections. The final values of the parameter μ are 6.5 × 10³, 1.2 × 10⁴ and 2.0 × 10⁴ for the subset numbers 1, 2, and 5, respectively. The reconstructed results are presented in Figs. 8–10, which show a satisfied stability of the imaging performance when the number of the OS is lower (OS = 1,2). However, the SSIM of OS-SB-5 is less stable than that of the noiseless projections and the corresponding MSE also becomes vibrating. The main reason is that there is no sufficient information in each subset especially for noise projections case. For the selection of the subset number, it is well known that the subsets should be balanced to make the pixel activity contributing equally to any subset. In addition, it is also imperative to ensure that each subset has sufficient statistical information [33]. We recommend that each subset should include 20-40 projection views for a better performance.

From the aforementioned experiments, it is noticed that all the reconstructed images from TV minimization algorithms are in a high precision inside the ROI, and this is consistent with the CS-based interior tomography theory. It can also be observed that when the subset number used for the SFT-OS-SART is small, the SB methods converge faster than the SFT-OS-SART. The use of the OS technique can further speedup the convergence of the SB methods. However, the increase of the subset number can somewhat decrease the image quality of the reconstructed images especially for the projections include noise.

4.3 Clinical application

To demonstrate the practical value of the proposed algorithm, cardiac images are reconstructed from a clinical projection data collected on a GE Discovery CT750 HD scanner at Wake Forest University Health Science. The scanning trajectory is a typical helix with a radius of 538.5 mm. After appropriate preprocessing, a set of 64-slices fan-beam sinograms are obtained. Over a 360-degree scanning range, 2200 views are uniformly acquired for each sinogram. For one projection, 888 detector elements are equi-angularly distributed, covering a field of view (FOV) of 249.2 mm in radius. In this experiment, we chose the 32-nd slice to evaluate the proposed algorithm. Figure 11 (a) gives the reconstructed image by the conventional filtered backprojection (FBP) method from all the 2200 views. To evaluate the OS-SB based interior reconstruction technique, the sinogram is manually truncated by discarding 300 detector elements on each side, covering a small ROI with a radius of 77.1 mm (see Fig. 11 (b)).

In our experiments, 1100 and 220 views are uniformly downsampled from all 2200 views, respectively. In the case of 1100 views, the images are reconstructed by the OS-SB method with 1, 2 and 5 subsets. For a fair comparison, the subset number 5 is selected for the OS-SART and STF-OS-SART. However, the OS-SART converges slowly when a small subset number is employed, and there is almost no difference between the OS-SART and STF-OS-SART algorithms. Therefore, a subset number 44 is also used to significantly accelerate the OS-SART and STF-OS-SART. The values of the parameters λ and μ of the OS-SB are summarized as in Table 2 for the clinical application. The reconstructed results after 20 iterations are shown in Fig. 12 in a display window of [–1000 HU, 1500HU]. While the first and third row’s images are reconstructed by different algorithms, the images on the second and fourth rows are the corresponding differences with respect to the reference image in Fig. 11 reconstructed by the FBP approach from all the projections. To quantitatively characterize the image quality, the MSE and SSIM were also calculated for the ROI inside the white circle as indicated in Fig. 11 (b).Figure 13 gives the MSE and SSIM curves with respect to the iteration numbers for different algorithms. Figure 14 shows the representative profiles along the white line indicated in Fig. 11 (b) (here we only list the best results of STF-OS-SART and OS-SB algorithms). From Figs. 11–14, we can see that the OS-SB-1 outperforms the OS-SART-5, and its performance is comparable to the STF-OS-SART-44. From Fig. 13, we can see that the OS-SB with 2 and 5 subsets have smaller MSEs compared with the OS-SART-44 and STF-OS-SART-44. However, the SSIM curves in Fig. 13 (b) are not entirely consistent to Fig. 13 (a). Except the aforementioned reasons for numerical simulations, additional reason is as follow. The OS-SB divides the projection matrix into several sub-matrices. For each sub-matrix, the sub problem is also modeled as a constrained problem and solved by the SB method. These subproblems are alternatively optimized. The SD based updating process assumes that the gradients of all the sub-objective functions are approximately the same. This assumption fails for the OS-SB when the subset number is greater than one.

For the case of 220 views, the images are also reconstructed by the OS-SB with 1, 2 and 5 subsets, respectively. For a fair comparison, we select 5 and 22 subsets for the OS-SART, and 22 subsets for the STF-OS-SART. The reconstructed results after 20 iterations are shown in Figs. 15–17. From Figs. 15 and 16, we can see that the differences between the interior reconstructions and reference image become smaller and smaller from (a’) to (f’), especially for the OS-SB-2 and OS-SB-5. We also notice that the SSIMs for the images reconstructed from 220 views are more unstable than those from 1100 views due to the reduction of the projections.

From the above experiments, one can see that the performances of the SB and STF-OS-SART are still acceptable for truncated clinical data. However, the performance of OS-SB with subsets is not as good as the numerical simulations from noise-free projections. This is due to several reasons. First, this realistic clinical image does not satisfy the ideal piecewise constant imaging model. Second, the measured projections include noise. Third, the DGT of realistic image may not sparse enough to meet the RIP condition in the CS theory for a set of downsampled projections. As an alternative way, other possible sparse transforms (such as wavelet transform [34], curvelet transform [35] or learned dictionary expression [36]) have been studied by peers.

5 Conclusion

In this paper, we have developed an OS-SB iterative algorithm for CS-based interior reconstruction and evaluated its performance with both simulated projections and realistic clinical data sets. To demonstrate the merits of the proposed method for ROI reconstruction without/with few-view projections, the OS-SART and STF-OS-SART are also implemented and compared with. Our experimental results show the SB approach outperforms the STF-OS-SART with a small subset number for interior tomography, especially for few-view case. The OS technique can further accelerate the convergence rate of the SB with additional noise. Therefore, a small subset number can be selected to compromise the reconstructed image quality and the convergence speed. In the near future, the GPU–based acceleration [37] can be implemented to further speedup the proposed algorithm especially for cone-bam interior tomography.

Footnotes

Acknowledgments

This work was partially supported by grants from the National Science Foundation of China (No. 61171179), the Science Foundation of Shanxi Province (No. 2012021011-2), as well as the National Science Foundation of USA under a CBET CAREER award (No. 1540898).

References

, Yu

and Wang

, Exact interior reconstruction with Cone-Beam CT, Int J Biomed Imaging 2007 (2007). Article ID. 10693.

, Yu

, Wei

and Wang

, A general local reconstruction approach based on a truncated hilbert transform, Int J Biomed Imaging 2007 (2007). Article ID. 63634.

, Yu

and Wang

, Exact interior reconstruction from truncated limited-angle projection data, Int J Biomed Imaging 2008 (2008). Article ID. 427989.

, Ye

and Wang

, Interior reconstruction using the truncated Hilbert transform via singular value decomposition, J X-Ray Sci Technol 16(4) (2008), 243–251.

Candes

E.J.

, Romberg

and Tao

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

Donoho

D.L.

, Compressed sensing, IEEE Trans Inf Theory 52(4) (2006), 1289–1306.

and Wang

, Compressed sensing based interior tomography, Phys Med Biol 54(9) (2009), 2791–2805.

Yang

, Yu

, Jiang

and Wang

, High-order total variation minimization for interior tomography, Inverse Probl 26(3) (2010). Article ID 035013.

, Wang

, Hsieh

, Entrikin

D.W.

, Ellis

, Liu

and Carr

J.J.

, Compressive Sensing-Based interior tomography: Preliminary clinical application, J Comput Assist Tomogr 35(6) (2011), 762–764.

10.

, Chen

and Zhou

, CT image reconstruction from sparse projections using adaptive TpV regularization, Comput Math Methods Med 2015 (2015). Article ID 354869.

11.

Zhang

and Xing

, Sequentially reweighted TV minimization for CT metal artifact reduction, Med Phys 40(7) (2013). Article ID 071907.

12.

Yin

, Osher

, Goldfarb

and Darbon

, Bregman iterative algorithms for l1-minimization with applications to compressed sensing, SIAM J Imaging Sci 1(1) (2008), 143–168.

13.

and Wang

, A soft-threshold filtering approach for reconstruction from a limited number of projections, Phys Med Biol 55(13) (2010), 3905–3916.

14.

Goldstein

and Osher

, The split bregman method for l1-regularized problems, SIAM J Imaging Sci 2(2) (2009), 323–343.

15.

Chang

, Li

, Chen

, Xiao

, Zhang

and Wang

, A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction, J X-Ray Sci Technol 21(2) (2013), 161–176.

16.

Vandeghinste

, Goossens

, Van Holen

, Vanhove

, Pizurica

, Vandenberghe

and Staelens

, Iterative CT reconstruction using shearlet-based regularization, IEEE Trans Nucl Sci 60(5) (2013), 3305–3317.

17.

Ramani

and Fessler

J.A.

, A Splitting-Based iterative algorithm for accelerated statistical X-Ray CT reconstruction, IEEE Trans Med Imaging 31(3) (2012), 677–688.

18.

Yang

, Qi

, Xu

, Zhen

, Lu

and Zhou

, Cone beam CT image iterative reconstruction based on Split-Bregman method, Nan Fang Yi Ke Da Xue Xue Bao 34(6) (2014), 783–786.

19.

Hudson

H.M.

and Larkin

R.S.

, Accelerated image reconstruction using ordered subsets of projection data, IEEE Trans Med Imaging 13(4) (1994), 601–609.

20.

Erdogan

and Fessler

J.A.

, Ordered subsets algorithms for transmission tomography, Phys Med Biol 44(11) (1999), 2835–2851.

21.

Kole

J.S.

and Beekman

F.J.

, Evaluation of the ordered subset convex algorithm for cone-beam CT, Phys Med Biol 50(4) (2005), 613–623.

22.

Zhao

, Yang

and Jiang

, A fast algorithm for high order total variation minimization based interior tomography, J X-Ray Sci Technol 23(3) (2015), 349–364.

23.

and Wang

, Finite detector based projection model for super resolution CT, J X-Ray Sci Technol 20(2) (2012), 229–238.

24.

Jiang

and Wang

, Convergence of the simultaneous algebraic reconstruction technique (SART), IEEE Trans Image Process 12(8) (2003), 957–961.

25.

Gordon

, Bender

and Herman

G.T.

, Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography, J Theor Biol 29(3) (1970), 471–481.

26.

Maximum likelihood reconstruction for emission tomography, IEEE Trans Med Imaging 1(2) (1982).

27.

Daubechies

, Defrise

and De Mol

, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Commun Pure Appl Math 57(11) (2004), 1413–1457.

28.

Wang

, Yin

and Zhang

, A Fast Algorithm for Image Deblurring with Total Variation Regularization. CAAM Technical Report TR07-10, 2007.

29.

Boyd

, Parikh

, Chu

, Peleato

and Eckstein

, Distributed optimization and statistical learning via the alternating direction method of multipliers, Found Trends Mach Learn 3(1) (2011), 1–122.

30.

Nocedal

Jorge

and Wright

Stephen J.

“Line Search Methods,” in Numerical Optimization, Springer New York, 2006, pp. 30–65.

31.

Wang

and Jiang

, Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART), J X-Ray Sci Technol 12(3) (2004), 169–177.

32.

Wang

, Bovik

A.C.

, Sheikh

H.R.

and Simoncelli

E.P.

, Image quality assessment: From error visibility to structural similarity, IEEE Trans Image Process 13(4) (2004), 600–612.

33.

Kadrmas

Dan J.

, Statistically regulated and adaptive EM reconstruction for emission computed tomography, IEEE Trans Nucl Sci 48(3) (2001), 790–798.

34.

Averbuch

, Lazar

and Israeli

, Image compression using wavelet transform and multiresolution decomposition, IEEE Trans Image Process 5(1) (1996), 4–15.

35.

Starck

J.-L.

, Candes

E.J.

and Donoho

D.L.

, The curvelet transform for image denoising, IEEE Trans Image Process 11(6) (2002), 670–684.

36.

, Yu

, Mou

, Zhang

, Hsieh

and Wang

, Low-Dose X-ray CT reconstruction via dictionary learning, IEEE Trans Med Imaging 31(9) (2012), 1682–1697.

37.

Liu

, Luo

and Yu

, GPU-based acceleration for interior tomography, IEEE Access 2 (2014), 757–770.