A fast iteration approach to undersampled cone-beam CT reconstruction

Abstract

BACKGROUND:

Due to large dimensional matrix multiplications, the existing iterative algorithms for cone beam computed tomography (CBCT) reconstruction often face problems of heavy computational workload and large volume of memory usage.

OBJECTIVE:

This study proposes and tests an iterative algorithm of 3DA-TVAL3 for fast reconstruction of CBCT images using undersampled measurement data and the reduced amount of computer memories.

METHODS:

In order to reduce computational workload and computer memories based on the sparsity of the CBCT measurement matrix, the proposed iterative algorithm applies elementwise scalar multiplications in the iterative computation to search for optimal solution. Through a number of tests on three different CT data sets with different number of projections, the reconstruction performance of the proposed algorithm is compared with that of two accelerated iterative algorithms and the conventional FDK algorithm.

RESULTS:

The visual and quantitative evaluations using the normalized mean square error, peak signal to noise ratio and structural similarity metrics demonstrated the faster computational time and the higher image quality of using the proposed 3DA-TVAL3 algorithm than using other conventional algorithms under comparison.

CONCLUSIONS:

The proposed 3DA-TVAL3 algorithm can perform efficient and fast computation of CBCT reconstruction using the reduced amount of computer memories.

Keywords

Cone-beam computed tomography optimization CT reconstruction

1 Introduction

Computed tomography (CT) is one of most significant innovations medical imaging in recent decades, which can provide fast and effective analysis and identification of lesions by providing 3D visual inspection of tissues and organs of patients. With the development of CT technologies, the on-board cone-beam CT (CBCT) attracts growing attention in radiation imaging and therapy for its high spatial resolution, easy operation and relatively low acquisition expense [1 –4].

Despite the high on-board volumetric imaging quality, there has been a concern of the excessive X-ray radiation does into the patient body, which is much higher than that of conventional diagnostic X-rays [5, 6]. It is expected that the CBCT scanning and imaging process can reduce the X-ray radiation does as much as possible while maintaining the quality of CT images. To reduce the number of X-ray projections in the CBCT scan directly reduces the radiation dose, which results in under sampled measurement data. In the literature, a type of filtered back-projection (FBP) algorithms, originally proposed by Feldkamp, Davis, and Kress (FDK) [7] and its derivatives [8, 9], have been widely applied to CBCT reconstruction. It is found that when the FDK type algorithms are applied to undersampled projections, the quality of reconstructed images are degraded dramatically. This is due to the incomplete information in the Fourier domain.

In contrast, another type of iterative reconstruction algorithms can provide better imaging quality using undersampled measurement data [10 –14]. This type of algorithms is to solve the linear matrix equation of the CT measurement data based on computational optimization of a specified cost function. Because the measurement data and the measurement matrix are of large dimension, the iteration algorithms carry out repeated large dimensional matrix multiplication operations for computation of the optimal reconstruction solution. These matrix multiplication operations take a large amount of computational workload and require a large amount of computer memories for storage of the matrix data.

This paper extends the result of [15] on computing the X-ray transform and adjoint matrices to present an algorithm for fast iterative CBCT reconstruction to deal with the large dimensional matrix multiplication problem. Based on that the rows and columns of the CBCT measurement matrix are sparse, the proposed algorithm performs direct elementwise operations on the measurement matrix multiplication to reduce the computational workload and usage of computer memories. Such elementwise computation does not require large dimensional measurement matrix multiplications and large amount of memory storages of the large dimensional attenuation coefficients. Only a small number of scalar multiplications and a small set of parameters of the structural and geometric information of the CBCT scan are used in each elementwise computation. Using the direct elementwise operation for measurement matrix multiplication, the 3D accelerated total variation minimization by augmented Lagrangian and alternating direction algorithm (3DA-TVAL3) is proposed for CBCT image reconstruction. The fast iterative 3DA-TVAL3 algorithm is proposed for circular CBCT reconstruction and is readily extendible to spiral CBCT reconstruction.

The rest of this paper is organized as follows. Section 2 on materials and methods presents formulations of the elementwise matrix multiplications and the augmented Lagrangian and alternating direction algorithm. Section 3 presents the reconstruction results of the proposed algorithm on three difference CT data sets at different numbers of projections in comparison with other accelerated algorithms and the conventional FDK algorithm. Finally, Section 4 is on discussion and conclusion.

2 Materials and methods

2.1 CT systematic description and acceleration

A circular CBCT scanning system basically consists of an X-ray source, the scanned object and a flat panel detector. The geometric configuration of scanning system in a Cartesian coordinate system is shown in Fig. 1. Suppose that the attenuation coefficient of reconstructed image is discretized into a 3D matrix with N_u = N_x × N_y × N_z voxels in the Cartesian x-y-z coordinate system, the size of each voxel is d_x × d_y × d_z, N_p represents the number of projections and N_a × N_b represents the size of each projection. Let u ∈ R^{N
_u} denote the vectorized image which can be written as

Fig. 1.

Configuration of the cone beam CT imaging geometry.

$u = {[u_{1}, u_{2}, . . ., u_{N_{u}}]}^{T} .$ (1)

As a result, the attenuation coefficient of a voxel indexed by (x, y, z) in the 3D object is the kth entry of the vector u, where $k = (z - 1) N_{x} N_{y} + (x - 1) N_{y} + y .$ (2)

Then the CBCT system can be formulated into the following linear equation $v = Au,$ (3)

where v ∈ R^{N
_u} is an N_v = N_p × N_a × N_b dimensional vector representing the projection measurements and A ∈ R^N_v×N_u is the measurement matrix.

Though iterative reconstruction methods often offer better image reconstruction quality than that of the conventional FBP type methods, there is a problem of heavy computational workload that the iterative methods have not been commonly applied to clinical diagnosis. The problem is caused by repetitive calculations of Au^t along with updates of the reconstructed image in each iteration step. To address this problem, this paper extends the result of [10] on computation of the X-ray transform and its adjoint to the iterative optimization for CBCT reconstruction. For reducing the computational workload and computer memories, our proposed iterative optimization algorithm carries out elementwise computation for updatingAu^t and its adjoint A^Tv^t, where u^t and v^t are the updated vectors at tth iteration step.

2.2 Direct computation of Au^t

Let A_ri ∈ R^1×N_vandA_cj ∈ R^{N
_u}, i = 1 … N_v, j = 1 … N_u, denote the row and column vectors of A, such that $A = {[A_{r 1}^{T} … A_{r N_{v}}^{T}]}^{T} = [A_{c_{1}} … A_{r N_{u}}] .$ Without loss of generality, let N_x = N_y, d_x = d_y = 1 and the isocenter be at the origin. Assume that a radiation ray intersects the whole 3D domain with start point (x^s, y^s, z^s) and end point (x^e, y^e, z^e). For construction of each A_ri of A, an observation is that the start and end points of a radiation ray, with respect to the scanned object, satisfy one and only one of the following conditions: |x^e - x^s| > |y^e - y^s| and |x^e - x^s| > |z^e - z^s|; |y^e - y^s| > |x^e - x^s| and |y^e - y^s| > |z^e - z^s|; or |z^e - z^s| > |x^e - x^s| and |z^e - z^s| > |y^e - y^s|.

For the simplicity and without loss of generality, consider the first case that the radiation ray satisfies |x^e - x^s| > |y^e - y^s| and |x^e - x^s| > |z^e - z^s|. The 3D image can be discretized into planes perpendicular to the x-axis in the x-y-z coordinate system, which are called y-z-planes. The radiation ray nontrivially intersects $(X min - X max + 1)$ y-z-planes with their index X satisfying $X min \leq X \leq X max$ , where $X min and X max$ are the greatest integers below min(x^s, x^e) and max(x^s, x^e), respectively. It may intersect 1, 2, or 3 voxels of each y-z-plane indexed by X with $X min \leq X \leq X max$ . Associated with the radiation ray, there is a row A_ri of A with entries representing intersection lengths of the radiation ray passing through the corresponding voxels of the $(X min - X max + 1)$ y-z-planes. For computing the intersection lengths of the voxels in each y-z-plane indexed by x, introduce notations $Y -, Y +, Z - and Z +$ as the greatest integers that are smaller than y^-, y⁺, z^- and z⁺, respectively, $k_{y} = \frac{y^{e} - y^{s}}{x^{e} - x^{s}}$ , $k_{z} = \frac{z^{e} - z^{s}}{x^{e} - x^{s}}$ , y^- = k_y (X - x^s) + y^s, y⁺ = k_y (X + 1 - x^s) + y^s, z^- = k_z (X - x^s) + z^s,z⁺ = k_z (X + 1 - x^s) + z^s, $r_{y} = \frac{max (Y^{-}, Y^{+}) - y^{-}}{y^{+} - y^{-}} and r_{z} = \frac{max (Z^{-}, Z^{+}) - z^{-}}{z^{+} - z^{-}} .$ _. y^-, y⁺, z^-andz⁺ are the y-coordinates and z-coordinates of the intersections of the beam and two boundaries in y-z-plane. r_y and r_z are the intersecting ratios of the beam along the y-direction and z-direction.

It follows from [15] that, under the conditions |x^e - x^s| > |y^e - y^s| and |x^e - x^s| > |z^e - z^s|, the intersection lengths of the voxels, denoted by l_(x,y,z) in each y-z-plane, indexed by x, can be computed by the formulae summarized in Table 1 for five different cases. It is shown that the intersection lengths l_(x,y,z) and the corresponding coordinates of the voxels in each y-z-plane can be configured based on the coordinate information of the start and end points, i.e. (x^s, y^s, z^s) and (x^e, y^e, z^e).

Table 1
Intersection coordinates and lengths of voxels in a y-z-plane

Case Number of intersection voxels Coordinate of intersection voxels Length of intersection Component $v_{i}^{x}$ of v_i contributed by the radiation ray intersecting with the voxels in the y-z-plane

$\begin{matrix} Y - = Y + \\ Z - = Z + \end{matrix}$ 1 (X, Y^-, Z^-) $l_{(X, Y^{-}, Z^{-})} = \sqrt{1 {+ k}_{y}^{2} + k z 2}$ $v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t}$

$\begin{matrix} Y - = Y + \\ Z - \neq Z + \end{matrix}$ 2 (X, Y^-, Z^-) $l_{(X, Y^{-}, Z^{-})} = r_{z} \sqrt{1 {+ k}_{y}^{2} + k z 2}$ $\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{-}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \end{matrix}$

(X, Y^-, Z⁺) $l_{(X, Y^{-}, Z^{+})} = (1 - r_{z}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

$\begin{matrix} Y - \neq Y + \\ Z - = Z + \end{matrix}$ 2 (X, Y^-, Z^-) $l_{(X, Y^{-}, Z^{-})} = r_{y} \sqrt{1 {+ k}_{y}^{2} + k z 2}$ $\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u^{t} (Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-} t \\ + l_{(X, Y^{+}, Z^{-})} u^{t} (Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+} \end{matrix}$

(X, Y⁺, Z^-) $l_{(X, Y^{+}, Z^{-})} = (1 - r_{y}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

$\begin{matrix} Y - \neq Y + \\ Z - \neq Z + \\ r_{y} > r_{z} \end{matrix}$ 3 (X, Y^-, Z^-) $l_{(X, Y^{-}, Z^{-})} = r_{z} \sqrt{1 {+ k}_{y}^{2} + k z 2}$ $\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{-}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{+}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+}}^{t} \end{matrix}$

(X, Y^-, Z⁺) $l_{(X, Y^{-}, Z^{+})} = (r_{y} - r_{z}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

(X, Y⁺, Z⁺) $l_{(X, Y^{+}, Z^{+})} = (1 - r_{y}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

$\begin{matrix} Y - \neq Y + \\ Z - \neq Z + \\ r_{y} \leq r_{z} \end{matrix}$ 3 (X, Y^-, Z^-) $l_{(X, Y^{-}, Z^{-})} = r_{y} \sqrt{1 {+ k}_{y}^{2} + k z 2}$ $\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{+}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+}}^{t} \\ + l_{(X, Y^{+}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+}}^{t} \end{matrix}$

(X, Y^-, Z⁺) $l_{(X, Y^{-}, Z^{+})} = (r_{z} - r_{y}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

(X, Y⁺, Z⁺) $l_{(X, Y^{+}, Z^{+})} = (1 - r_{z}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

Case	Number of intersection voxels	Coordinate of intersection voxels	Length of intersection	Component $v_{i}^{x}$ of v_i contributed by the radiation ray intersecting with the voxels in the y-z-plane
$\begin{matrix} Y - = Y + \\ Z - = Z + \end{matrix}$	1	(X, Y^-, Z^-)	$l_{(X, Y^{-}, Z^{-})} = \sqrt{1 {+ k}_{y}^{2} + k z 2}$	$v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t}$
$\begin{matrix} Y - = Y + \\ Z - \neq Z + \end{matrix}$	2	(X, Y^-, Z^-)	$l_{(X, Y^{-}, Z^{-})} = r_{z} \sqrt{1 {+ k}_{y}^{2} + k z 2}$	$\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{-}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \end{matrix}$
		(X, Y^-, Z⁺)	$l_{(X, Y^{-}, Z^{+})} = (1 - r_{z}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$
$\begin{matrix} Y - \neq Y + \\ Z - = Z + \end{matrix}$	2	(X, Y^-, Z^-)	$l_{(X, Y^{-}, Z^{-})} = r_{y} \sqrt{1 {+ k}_{y}^{2} + k z 2}$	$\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u^{t} (Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-} t \\ + l_{(X, Y^{+}, Z^{-})} u^{t} (Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+} \end{matrix}$
		(X, Y⁺, Z^-)	$l_{(X, Y^{+}, Z^{-})} = (1 - r_{y}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$
$\begin{matrix} Y - \neq Y + \\ Z - \neq Z + \\ r_{y} > r_{z} \end{matrix}$	3	(X, Y^-, Z^-)	$l_{(X, Y^{-}, Z^{-})} = r_{z} \sqrt{1 {+ k}_{y}^{2} + k z 2}$	$\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{-}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{+}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+}}^{t} \end{matrix}$
		(X, Y^-, Z⁺)	$l_{(X, Y^{-}, Z^{+})} = (r_{y} - r_{z}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$
		(X, Y⁺, Z⁺)	$l_{(X, Y^{+}, Z^{+})} = (1 - r_{y}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$
$\begin{matrix} Y - \neq Y + \\ Z - \neq Z + \\ r_{y} \leq r_{z} \end{matrix}$	3	(X, Y^-, Z^-)	$l_{(X, Y^{-}, Z^{-})} = r_{y} \sqrt{1 {+ k}_{y}^{2} + k z 2}$	$\begin{matrix} v_{i}^{x} = l_{(X, Y^{-}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{-}}^{t} \\ + l_{(X, Y^{+}, Z^{-})} u_{(Z^{-} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+}}^{t} \\ + l_{(X, Y^{+}, Z^{+})} u_{(Z^{+} - 1) N_{x} N_{y} + (X - 1) N_{y} + Y^{+}}^{t} \end{matrix}$
		(X, Y^-, Z⁺)	$l_{(X, Y^{-}, Z^{+})} = (r_{z} - r_{y}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$
		(X, Y⁺, Z⁺)	$l_{(X, Y^{+}, Z^{+})} = (1 - r_{z}) \sqrt{1 {+ k}_{y}^{2} + k z 2}$

Under the conditions |x^e - x^s| > |y^e - y^s| and |x^e - x^s| > |z^e - z^s| and with a given image update u^t, the intersection lengths obtained as in Table 1 can be extended to compute entries of the measurement vector v_i = A_riu^t, for 1 ≤ i ≤ N_v. For each y-z-plane, indexed by x, the radiation ray intersects at most three of its voxels. With respect to the five cases as listed in Table 1, the component $v_{i}^{x}$ of v_i contributed by the radiation ray intersecting with the voxels in the y-z-plane can be computed as summarized in the last column of Table 1. It follows that v_i = A_riu^t can be computed as $v_{i} = \sum_{X - \leq x \leq X^{+}} v_{i}^{x} .$ (4)

This result under the conditions |x^e - x^s| > |y^e - y^s| and |x^e - x^s| > |z^e - z^s| can be generalized as follows $v_{i} = {\begin{matrix} \sum_{X^{-} \leq x \leq X^{+}} v_{i}^{x}, if | x^{e} - x^{s} | > | y^{e} - y^{s} |, | x^{e} - x^{s} | > | z^{e} - z^{s} |, \\ \sum_{Y^{-} \leq y \leq Y^{+}} v_{i}^{y}, if | y^{e} - y^{s} | > | x^{e} - x^{s} |, | y^{e} - y^{s} | > | z^{e} - z^{s} |, 1 \leq i \leq N_{v}, \\ \sum_{Z^{-} \leq z \leq Z^{+}} v_{i}^{z}, if | z^{e} - z^{s} | > | x^{e} - x^{s} |, | z^{e} - z^{s} | > | y^{e} - y^{s} |, \end{matrix}$ (5)

where the definitions of the x-z-planes and x-y-planes follow from that of the y-z-planes and the computations of $v_{i}^{y}$ , for Y_min ≤ Y ≤ Y_max, and $v_{i}^{z}$ , for Z_min ≤ Z ≤ Z_max, follow from the similar computational steps for $v_{i}^{x}$ .

As a result, for a given image update u^t, equation (4), based on the technical results summarized in Table 1, provides an approach to the direct computation of Au^t. A special characteristic of this approach is that it only requires geometric information of the scanning radiation rays in terms of the start and end points (x^s, y^s, z^s) and (x^e, y^e, z^e) but does not involve explicit formation of the large dimensional matrix A and matrix multiplication.

2.3 Direct Computation of A^Tv^t

For elementwise computation of the adjoint transform A^Tv^t, it is noted that each adjoint element $u_{j}^{t} = A_{cj}^{T} v^{t}$ represents the sum of lengths of the radiation rays through the jth voxel of the scanned object, for 1 ≤ j ≤ N_u. Let (x_cj, y_cj, z_cj) be the coordinate of the jth voxel, d_a × d_b the size of each detector sensor, (n_a, n_b) the 2D index of the detector panel of size N_a × N_b and n_p × dφ the projection angle of the n_pth projection with n_p ∈ [0, N_p - 1] and dφ being the angle variance between two projections. The distance from the source to the isocenter denoted by S_O and that from the source to the center of detection panel is denoted by S_D.

For a fixed scan projection on a voxel of the object, there exist radiation beams passing through the eight corners of the voxel as shown in Fig. 2. The straight lines connecting the intersection points of these beams on the detector panel form a closed region on the detector panel. And the radiation beams through the voxel are only projected onto the detectors within this closed region. At each projection angle n_p × dφ, the coordinate of the voxel at (x_cj, y_cj, z_cj) is transformed to a new rotated coordinate system (x_rj, y_rj, z_cj), so the detector panel is perpendicular to the x_r-axis of the transformed coordinate system. As a result $\begin{matrix} x_{r} = x_{c} cos (n_{p} \times d φ) + y_{c} sin (n_{p} \times d φ) & and & y_{r} = - x_{c} sin (n_{p} \times d φ) + y_{c} cos (n_{p} \times d φ) \end{matrix} .$ (6)

Fig. 2.

The projected X-ray beams through a voxel onto a region of the detector panel. S is the source of x-ray and D is the detector panel. The inner cube represents the jth voxel. The ball encloses the voxel and the outer cube encloses the ball. (a) the 3D view, (b) the 2D view in x-y plane.

To determine the enclosed region on the detector panel, a ball that encloses the voxel is constructed and an upright cube is formed to enclose the ball as shown in Fig. 2. Then boundaries on the detector panel, in terms of the panel detector indices can be found as follows.

$\begin{matrix} n_{aj, min} = \frac{(y_{rj} - 1) S_{D}}{(x_{rj} + S_{O}) d_{a}} + \frac{N_{a}}{2}, n_{aj, max} = \frac{(y_{rj} + 1) S_{D}}{(x_{rj} + S_{O}) d_{a}} + \frac{N_{a}}{2}, \\ n_{bj, min} = \frac{(z_{cj} - Δ) S_{D}}{(x_{rj} + S_{O}) d_{b}} + \frac{N_{b}}{2}, n_{bj, max} = \frac{(z_{cj} + Δ) S_{D}}{(x_{rj} + S_{O}) d_{b}} + \frac{N_{b}}{2}, Δ = \sqrt{1 + \frac{{d_{z}}^{2}}{2}} \end{matrix}$ (7)

In fact, for each voxel of the object, there are only a small number of detectors that collect the radiation beams through the voxel. The intersection length of each of these beams through the fixed voxel can be computed according to the start and end points of each of these beams and following the formulae as given in Table 1. As a result, the elementwise computation of $u_{j}^{t} = A_{cj}^{T} v^{t}$ , j = (z_cj - 1) N_xN_y + (x_rj - 1) N_y + y_rj, is $u_{j}^{t} = \sum_{\begin{matrix} n_{aj, min} \leq n_{a} \leq n_{aj, max} \\ \begin{matrix} n_{bj, min} \leq n_{b} \leq n_{bj, max} \\ 0 \leq n_{p} \leq N_{p} - 1 \end{matrix} \end{matrix}} l_{(x_{rj}, y_{rj}, z_{cj})}^{(n_{p}, n_{a}, n_{b})} v_{(n_{p} n_{b} n_{a} + n_{b} N_{a} + n_{a})}^{t},$ (8)

where $l_{(x_{rj}, y_{rj}, z_{cj})}^{(n_{p}, n_{a}, n_{b})}$ denotes the intersection length through the voxel at (x_cj, y_cj, z_cj) and projected onto the detector indexed by (n_a, n_b).

2.4 The iteration algorithm 3DA-TVAL3

The subsection proposes a fast iteration algorithm for CBCT reconstruction using the elementwise computation of elementwise computation of Au^t and A^Tv^t in the iterative updates. It is known that CT images have sparse gradient-magnitude images (GMI) and minimizing the image total variation (TV) norm, which is the l₁-norm of the GMI, subject to data constraints can produce satisfactory results from sparse data samples [11]. Let D_iu ∈ R³ denote the discrete gradient of the ith voxel of u. Given the measurement data vector v satisfying the linear matrix equation (3), the TV regularization model for the cone-beam CT reconstruction can be formulated into the following constrained optimization problem:

$\begin{matrix} min_{u} & TV (u) = \sum_{i} {∥ D_{i} u ∥}_{2} & s . t . & v = Au \end{matrix} .$ (9)

To separate the non-differentiable l₁-norm term to reduce the computational workload in the optimization, g_i ∈ R³ is introduced as an equivalent of discrete gradient to split variables. Then the TV regularization model is reformulated as $\begin{matrix} min_{g_{i}, x} & \sum_{i} {∥ g_{i} ∥}_{2} & s . t . & v = Au, D_{i} u = g_{i} \end{matrix}, \forall i .$ (10)

For reconstruction from noise-corrupted measurements, the constrained optimization problem is formulated into an unconstrained optimization problem by introducing the following augmented Lagrangian function $min_{g_{i}, u} L_{A} (u, g_{i}, f_{i}, λ) = \sum_{i} ({∥ g_{i}, ∥}_{2} - f_{i}^{T} (D_{i} u - g_{i}) + \frac{α_{i}}{2} {∥ D_{i} u - g_{i} ∥}_{2}^{2}) - λ^{T} (Au - v) + \frac{τ}{2} {∥ Au - v ∥}_{2}^{2},$ (11)

where f_i ∈ R³, λ ∈ R^{N
_v} are the Lagrangian multipliers, τ is the compensation factor and α_i > 0 is the penalty parameter.

To compute the solution for the optimal minimization of the augmented Lagrangian function (11), the Alternating Direction Method (ADM) [16] is applied here to perform the following update of $(g_{i}^{t}, u^{t}, {f_{i}}^{t}, λ^{t})$ at each iteration step t. ${\begin{matrix} {g_{i}}^{t + 1} = Shrink 2 (D_{i} u^{t} - \frac{{f_{i}}^{t}}{α_{i}}, \frac{1}{α_{i}}), \forall i \\ u^{t + 1} = u^{t} - ɛ t d^{t} (u^{t}) \\ {f_{i}}^{t + 1} = {f_{i}}^{t} - η α_{i} (D_{i} u^{t + 1} - {g_{i}}^{t + 1}), \forall i \\ λ^{t + 1} = λ^{t} - η τ (A u^{t + 1} - v) \end{matrix}$ (12)

where η is a step length in the update of $λ$ , and it can guarantee the convergence of ADM when $η \in (0, \frac{\sqrt{5} + 1}{2})$ verified by theory [17, 18]. g_i is updated by the shrinkage formula which is defined as ${g_{i}}^{t + 1} = Shrink 2 (D_{i} u^{t} - \frac{{f_{i}}^{t}}{α_{i}}, \frac{1}{α_{i}}) = max {{∥ D_{i} u^{t} - \frac{{f_{i}}^{t}}{α_{i}} ∥}_{2} - \frac{1}{α_{i}}, 0} \cdot \frac{D_{i} u^{t} - \frac{{f_{i}}^{t}}{α_{i}}}{{∥ D_{i} u^{t} - \frac{{f_{i}}^{t}}{α_{i}} ∥}_{2}} .$ (13)

u is updated by the one-step steepest descent method. With fixed g_i^t+1, f_i^tand $λ$ ^t, L_A (u, g_i, f_i, λ) is quadratic with respect to u and its gradient can be easily derived as $d^{t} (u^{t}) = \sum_{i} [α_{i} {D_{i}}^{T} (D_{i} u^{t} - {g_{i}}^{t + 1}) - {D_{i}}^{T} {f_{i}}^{t}] + τ A^{T} (A u^{t} - v) - A^{T} λ^{t} \begin{matrix} . \end{matrix}$ (14)

Based on the above formulated optimization problem, the computational procedure of the proposed 3DA-TVAL3 is as follows.

Initialization: Input N_x, N_y, N_z, d_x, d_y, d_z and geometry parameters N_p, N_a, N_b, d_a, d_b, S_D, S_O of the CT model.

Input projections v and initialize output u, the Lagrangian multipliers f_i, $λ$ , iteration coefficient α_i, τ and stopping criterion tol, tol _ inn and maximum iteration number T.

t : =1

While t ≤ T Do

if $\frac{∥ u^{t} - u^{t - 1} ∥}{∥ u^{t - 1} ∥} > tol_inn$

calculate g^t+1 using Equation (13)

calculate step length d^t (u^t) using Equation (14)

update u^t+1 by one-step steepest descent method: u^t+1 = u^t - ɛ^td^t (u^t)

else if $\frac{∥ u^{t} - u^{t - 1} ∥}{∥ u^{t - 1} ∥} > tol$

update Lagrangian multipliers f_i, $λ$ using Equation (12)

calculate current augmented Lagrangian function using Equation (11)

else

output u

end if

t = t + 1

end while

2.5 Quantitative performance evaluation

For quantitative analysis, three error metrics were used. The first is the normalized mean square error (NMSE) defined below which is an index of error between the reference image u₀ and the reconstructed image u. $NMSE = \frac{{∥ u - u_{0} ∥}_{2}^{2}}{{∥ u 0 ∥}_{2}^{2}} .$ (15)

The second one is the peak signal to noise ratio (PSNR) which can measure the accuracy of reconstruction subject to noise. The PSNR between the reference image u₀ and the reconstructed image u is defined as: $PSNR = 10 {log}_{10} (max (u) / MSE (u_{0}, u))$ (16)

where max(u) is the maximum value of reference image and the mean square error (MSE) is defined as: $MSE (u_{0}, u) = \frac{1}{N_{x} N_{y} N_{z}} \sum_{i = 1}^{N_{x}} \sum_{j = 1}^{N_{y}} \sum_{k = 1}^{N_{z}} {[u_{0} (i, j, k) - u (i, j, k)]}^{2}$ (17)

The PSNR value approaches infinity as the MSE approaches zero. A higher PSNR value indicates higher reconstruction accuracy and, at the other side of the scale, a smaller value of the PSNR implies higher reconstruction error.

The third one is the structural similarity index (SSIM) which is a measure of similarity between the reference image u₀ and reconstructed image u in terms of structure and visual perception and is defined as $SSIM (u, u_{0}) = \frac{2 \bar{u} \cdot \bar{u_{0}} (2 σ_{u u_{0}} + c_{2})}{({\bar{u}}^{2} + {\bar{u_{0}}}^{2} + c_{1}) (σ_{u}^{2} + σ_{u_{0}}^{2} + c_{2})}$ (18)

where $\bar{u}$ and $\bar{u_{0}}$ represent mean values of u and u₀, respectively, $σ_{u}^{2}$ and $σ_{u_{0}}^{2}$ represent variances of u and u₀, respectively, σ_{uu
₀} represents covariance of u and u₀ and c₁ and c₂ are constant coefficient which can be set according to [19]. A higher SSIM value in general presents better reconstruction performance.”

3 Results

3.1 Experiment setup

The proposed 3DA-TVAL3 for CBCT reconstruction was applied, respectively, to a 3D digital Sheep-Logan phantom, a set of real CT data downloaded from the cancer imaging archive (TCIA) [20] and a set of scanned Catphan 500 physical phantom data acquired by a CBCT system of our institute. In the first two cases, the size of the Sheep-Logan phantom is 128×128×128 voxels, that of the real CT is 512×512×138 voxels, the source to detector distance is 1000 mm, that to the rotation center is 550 mm, the size of detection panel is 512×512 mm² and the size of each projection is 1024×1024 voxels. Since the Sheep-Logan phantom is smaller than the real CT data, only a quarter of the panel detector for faster reconstruction was used for simulation of the Sheep-Logan phantom. In the third case of the Catphan 500 physical phantom scan, the source to detector distance is 952 mm, that to the rotation center is 593 mm, the size of detection panel is 430×71.73 mm², the size of each projection is 1440×240 voxels, the number of projections is 502 with uniformly spaced angle increment. All the tests were implemented in the Matlab R2014a programming environment. The codes were implemented on a Win7 system using Intel(R) Xeon(R) CPU E5-2650 0 @ 2.00 GHz and 48.00 GB memory. The elementwise computation of transform and adjoint transform in section 2.2 and 2.3 were accelerated by a NIVIDIA Tesla K20c GPU.

3.2 Reconstruction of digital Sheep-Logan phantom

The 3D Sheep-Logan phantom was first tested, respectively, from 360, 180, 90, 60 projections along a circular orbit uniformly distributed in [0, 2π]. The reconstruction results of the proposed 3DA-TVAL3 were compared with that by FDK method. The mid slices of original phantom and reconstructed images in the x-y plane are shown in Fig. 3. With the decrease of projections, the performance of reconstruction results using the proposed 3DA-TVAL3 method can well preserve the critical information while the results by FDK method display considerable noise and artifacts.

Fig. 3.

The comparison of mid slice of 3D Sheep-Logan phantom reconstructed by FDK and 3DA-TVAL3. (a) the origin phantom; (b)∼(e) the FDK reconstruction results using 360, 180, 90, 60 projections respectively in the left column; (f)∼(i) the 3DA-TVAL3 reconstruction results using 360, 180, 90, 60 projections respectively in the right column.

To further illustrate the edge information, Fig. 4 shows 1D profiles of the original phantom and reconstructed images in x and y directions. It shows that results by 3DA-TVAL3 can preserve the image fidelity with highly undersampled projections, which outperform that of the FDK method.

Fig. 4.

1D profiles of the Sheep-Logan phantom in x and y directions which reconstructed by 3DA-TVAL3 and FDK. (a), (b) using 360 projections; (c), (d) using 180 projections; (e), (f) using 90 projections; (g), (h) using 60 projections.

To test the denoising effect of the proposed 3DA-TVAL3, additive Gaussian random noise matrices of zero-mean and variances 0.01, 0.05 and 0.10 were applied, respectively, to the projection data set acquired from 90 views. For each case of different variance, the noise matrix is produced using the normally distributed random number generation function of MATLAB. The reconstruction is repeatedly tested for 5 times with test using a newly generated noise matrix. As the ASD-POCS presented in [11] is a popularly known TV regularization method and EM-TV [21 –23] is a recently popular TV regularization method, their performance was compared with our proposed method. The averaged NMSEs, PSNRs and SSIMs of 3DA-TVAL3, ASD-POCS and EM-TV reconstructions subject to different external noise data are shown in Table 2. It is shown that the averaged NMSE is increased and that of PSNR and SSIM are reduced as the variance value of noise increased. All the averaged PSNRs of the noisy images reconstructed by 3DA-TVAL3 are higher than 40, which indicates that the differences between the reconstructed images and reference cannot be distinguished by human eye. In contrast, the averaged PSNRs of other TV methods are considerably smaller. So the proposed 3DA-TVAL3 algorithm outperforms the other algorithms in noise suppression.

Table 2

Averaged NMSEs, PSNRs and SSIMs of 3DA-TVAL3, ASD-POCS and EM-TV reconstructions with different external noise levels in the case of 90 projections

Noise Level	Method	NMSE	PSNR	SSIM
0	ASD-POCS	0.1512	30.3424	0.970944
	EM-TV	0.0675	37.2600	0.973417
	3DA-TVAL3	0.0070	57.1100	0.999636
1%	ASD-POCS	0.1499	30.3278	0.967859
	EM-TV	0.0694	37.0218	0.969777
	3DA-TVAL3	0.0075	56.3343	0.997905
5%	ASD-POCS	0.1514	30.2422	0.959136
	EM-TV	0.0718	36.7201	0.960550
	3DA-TVAL3	0.0107	53.2594	0.991513
10%	ASD-POCS	0.1528	30.1616	0.949359
	EM-TV	0.0738	36.4863	0.950665
	3DA-TVAL3	0.0137	51.1374	0.983745

The ASD-POCS and EM-TV algorithms were also accelerated by the proposed computation of the transform and adjoint transform in Subsections 2.2 and 2.3. The numbers of iterations of the ASD-POCS and EM-TV algorithms were set to be the same as that of the 3DA-TVAL3 to compare their convergence performance. Since the performance of reconstructed images degrades with reduction of number of projections, the reconstruction performance of these algorithms in the worst case of 60 projections is evaluated. The NMSE, PSNR and SSIM values of the three algorithms are shown in Table 3. It is shown that the 3DA-TVAL3 algorithm has the lowest NMSE and the highest PSNR and SSIM values with about the same reconstruction time as that of the other algorithms.

Table 3

Computational time, NMSEs, PSNRs and SSIMs of 3DA-TVAL3, ASD-POCS and EM-TV reconstructions with fixed number of iterations

Method	No. of iterations	Computational time	NMSE	PSNR	SSIM
ASD-POCS	165	663s	0.1245	27.5110	0.974100
EM-TV	165	613s	0.0351	42.9449	0.987964
3DA-TVAL3	165	642s	0.0093	54.4807	0.999544

Another experiment is carried out to further evaluate the reconstruction speed of the proposed 3DA-TVAL3 algorithm in comparison with that of the ASD-POCS and EM-TV algorithms. The three algorithms were set to terminate at the same NMSE value of 0.0134 and their numbers of iterations and computation times are recorded and shown in Table 4. It can be seen that, compared with ASD-POCS and EM-TV, the 3DA-TVAL3 algorithm has the least number of iterations and computational time.

Table 4

Number of iterations and computational time of 3DA-TVAL3, ASD-POCS and EM-TV reconstructions with fixed terminating NMSE

Method	No. of iterations	Computational time	NMSE	PSNR	SSIM
ASD-POCS	877	3583s	0.0134	51.3105	0.998416
EM-TV	952	3427s	0.0134	51.2567	0.994089
3DA-TVAL3	152	522s	0.0134	51.3206	0.998700

It is important to mention the NMSE value of the EM-TV algorithm eventually converged to NMSE = 0.0134 in this test. It changed little even if the algorithm further iterated to more than 952 times. Whereas the ASD-POCS and 3DA-TVAL3 algorithms could to converge to NMSE = 0.01 and even less. For this reason, their computational times were compared under the condition that their NMSE arrived at the same level 0.0134.

The mid slices of reconstructed images in transverse plane, sagittal plane and frontal plane are shown in Figs. 5–7, separately. Some boundary details are quite fuzzy in the reconstructed results of ASD-POC and EM-TV compared with that of the proposed 3DA-TVAL3. And there is little difference between the 3DA-TVAL3 reconstructed results and the references. Meanwhile, there are no significant artifacts caused by sparse sampling are observed in the boundary slices. All these can show that the proposed 3DA-TVAL3 algorithm can effectively suppress the noise and preserve the structural information in the reconstructed images.

Fig. 5.

Mid slices in transverse plane from (a) The origin phantom; (b) ASD-POCS reconstruction; (c) EM-TV reconstruction; (d) 3DA-TVAL3 reconstruction. The number of iteration is 165 and the number of projections is 60. All the images are displayed with the same window.

Fig. 6.

Mid slices in sagittal plane from (a) The origin phantom; (b) ASD-POCS reconstruction; (c) EM-TV reconstruction; (d) 3DA-TVAL3 reconstruction. The number of iteration is 165 and the number of projections is 60. All the images are displayed with the same window.

Fig. 7.

Mid slices in frontal plane from (a) The origin phantom; (b) ASD-POCS reconstruction; (c) EM-TV reconstruction; (d) 3DA-TVAL3 reconstruction. The number of iteration is 165 and the number of projections is 60. All the images are displayed with the same window.

Since the ASD-POCS and EM-TV algorithms have to use more computing time to carry out further iterations for possible improvement of the reconstruction quality, our proposed 3DA-TVAL3 algorithm certainly requires less computational time to produce the same reconstruction quality as that of the ASD-POCS and EM-TV algorithms. Hence our proposed 3DA-TVAL3 can perform faster reconstruction that performed by the accelerated ASD-POCS and EM-TV algorithms.

3.3 Reconstruction of downloaded CT data

The proposed 3DA-TVAL3 algorithm, together with the ASD-POCS and EM-TV algorithms, was applied to a set of real CT data downloaded from the cancer imaging archive for reconstruction. The quantitative evaluation results of the three algorithms, in terms of NMSE, PSNR and SSIM, are listed in Table 5. It is shown that the NMSE of 3DA-TVAL3 is significantly smaller than that of ASD-POCS and EM-TV all along for all tested projection numbers. The PSNR of 3DA-TVAL3 is about twice higher than that of ASD-POCS and EM-TV for all cases. With reduction of projection numbers, the SSIM of 3DA-TVAL3 remained almost unchanged while the SSIMs of ASD-POCS and EM-TV drop rapidly. These results can show that the 3DA-TVAL3 algorithm has advantages in noise suppression and structural information preservation in the case of undersampled projections. To illustrate visual superiority of the reconstruction by 3DA-TVAL3, the mid slices of the reconstructed images using ASD-POCS, EM-TV and 3DA-TVAL3 algorithms are shown in Fig. 8. It is shown that the boundary and detail of ASD-POCS reconstructed images are blurred. The images reconstructed by EM-TV preserve more structural information but still display considerable streaking artifacts as the projection number is reduced. In contrast, most structural information is preserved and few streaking artifacts exist in the images reconstructed by 3DA-TVAL3 with small projection numbers. These results can demonstrate the visual quality of the reconstruction by the proposed 3DA-TVAL3 algorithm subject to under sampled projections.

Fig. 8.

Comparison of mid slice of reconstructed images by the three algorithm. (a) the origin phantom; (b)∼(e) reconstructed images by ASD-POCS using 360, 180, 90, 60 projections respectively; (f)∼(i) reconstructed images by EM-TV using 360, 180, 90, 60 projections respectively; (j)∼(m) reconstructed images by 3DA-TVAL3 using 360, 180, 90, 60 projections respectively.

Table 5

NMSEs, PSNRs and SSIMs of 3DA-TVAL3, ASD-POCS and EM-TV reconstruction with different number of projections

No. Projections	Method	NMSE	PSNR	SSIM
360	ASD-POCS	0.1062	31.9949	0.920034
	EM-TV	0.0765	34.8487	0.939198
	3DA-TVAL3	0.0506	62.5174	0.999492
180	ASD-POCS	0.1071	31.9274	0.917322
	EM-TV	0.0786	34.6102	0.935334
	3DA-TVAL3	0.0567	61.5382	0.999338
90	ASD-POCS	0.1158	31.2427	0.897118
	EM-TV	0.0944	33.0217	0.918187
	3DA-TVAL3	0.0698	59.7272	0.999106
60	ASD-POCS	0.1273	30.4218	0.867041
	EM-TV	0.1066	31.9651	0.904005
	3DA-TVAL3	0.0790	58.6490	0.998712

3.4 Reconstruction of scanned CT data

The proposed 3DA-TVAL3 algorithm was also applied to the reconstruction of a set of Catphan 500 physical phantom data scanned by a CBCT system of our institute. The reconstruction results were compared with that the conventional FDK algorithm. Figure 7 shows the spatial resolution slices of the 3D Catphan 500 phantom reconstructed by the FDK and 3DA-TVAL3 algorithms, respectively. With the reduction of projection number to 1/7 (= 72 projections) of the full projection number, reasonable reconstruction quality was achieved by the proposed 3DA-TVAL3 as shown in Fig. 9 (b). In comparison, images reconstructed by the FDK methods using 502 projections and 72 projections are shown in Fig. 9 (a) and Fig. 9 (c), respectively, where the better reconstruction result of the proposed 3DA-TVAL3 algorithm can be clearly seen.

Fig. 9.

Spatial resolution slices of the reconstructed Catphan 500 physical phantom using (a) the FDK reconstruction from 502 projections; (b) the 3DA-TVAL3 reconstruction from 72 projections; (c) the FDK reconstruction from 72 projections.

4 Discussion and Conclusion

4.1 Discussion

A distinctive characteristic of the proposed 3DA-TVAL3 algorithm is that it performs elementwise computation of the matrix multiplication operations in each iterative updating. As formulated in Subsections 2.1–2.3, the elementwise computation in the iterations only uses a few geometrical parameters of the scan system without requiring storage of the large dimensional X-ray transform matrix. This directly leads to reduction of a significant amount of memories in the reconstruction processor. As a result, the RAM requirement of our proposed algorithm is about the same order of the reconstructed image size. For example, for reconstruction of the downloaded CT data of size 512×512×138 in Subsection 3.3, it was sufficient with 1GB of RAM in the GPU for implementing the 3DA-TVAL3 reconstruction.

The formulations of the elementwise operations in Subsections 2.1–2.2, as well as that in the Appendix, can indicate that the number of scalar multiplications is of the same order as the sparse number of voxel intersections. This implies that the proposed 3DA-TVAL3 algorithm carries out least number of scalar multiplications in each iteration.

Evaluations of the reconstruction speed of the proposed 3DA-TVAL3 algorithm were in comparison with the ASD-POCS and EM-TV algorithms. In the evaluations, the elementwise scalar multiplication operations on the transform and adjoint matrix multiplications were also implemented in the ASD-POCS and EM-TV algorithms. This is to make fair evaluations of the reconstruction speeds that the amounts of RAM usage of these algorithms were the same.

The fast reconstruction speed of the proposed 3DA-TVAL3 algorithm is under the condition that it is implemented with a small amount of RAM in the processor. There exist other reconstruction algorithms which may have faster reconstruction speed but may require more memory implementation in the processor. It is noted that not every reconstruction algorithm is feasible for direct implementation of the elementwise operations of the 3DA-TVAL3 algorithm. For example, the accelerated fast iterative shrinkage thresholding algorithms in [24] were based on preconditioning of the X-ray transform matrix. The proposed formulations for elementwise operations in Subsections 2.1–2.3 are not directly applicable to the preconditioned X-ray transform matrix.

4.2 Conclusion

An iterative 3DA-TVAL3 algorithm for CBCT reconstruction has been proposed in this paper. This algorithm is formulated as an augmented Lagrangian and alternating direction optimization problem and it performs elementwise computation of the matrix multiplication operations in each iterative update. A distinctive feature of the elementwise computation in the iterations is that it only uses geometrical parameters of the scan systems without storage of the X-ray transform matrix. This results in significant reduction of memories in the processor. Because of the sparse nature of the row and column vectors of the measurement matrix, each elementwise computation only involves a small sum of scalar multiplications, so it can result in efficient and fast computation of the iteration algorithm. Extensive qualitative and quantitative evaluations of the proposed 3DA-TVAL3 reconstruction algorithm were carried out in comparison with the conventional FDK and accelerated ASD-POCS and EM-TV algorithms. The evaluation results can demonstrate that the 3DA-TVAL3 reconstruction algorithm has advantages in memory reduction, noise suppression, reconstruction speed and accuracy.

Footnotes

Acknowledgments

This work was supported by the National Natural Science Foundation of China [81571772] and [61801475], by Jiangsu Planned Projects for Postdoctoral Research Funds [2018K180C] and by Key Science and Technology Plan of Jiangsu [BE2017671].

References

Yang ,

Schreibmann ,

Li ,

Wang and

Xing , Evaluation of on-board kV cone beam CT (CBCT)-based dose calculation, Physics in Medicine & Biology 52 (2007), 685–705.

Xing ,

Thorndyke ,

Schreibmann ,

Yang ,

T.F.

Li ,

G.Y.

Kim ,

Luxton and

Koong , Overview of image-guided radiation therapy, Medical Dosimetry 31 (2006), 91–112.

Lee ,

Q.T.

Le and

Xing , Retrospective IMRT dose reconstruction based on cone-beam CT and MLC log-file, International Journal of Radiation Oncology*Biology*Physics 70 (2008), 634–644.

D.A.

Jaffray ,

J.H.

Siewerdsen ,

J.W.

Wong and

A.A.

Martinez , Flat-panel cone-beam computed tomography for image-guided radiation therapy, International Journal of Radiation Oncology*Biology*Physics 53 (2002), 1337–1349.

M.J.

Murphy ,

Balter ,

J.A.

BenComo ,

I.J.

Das ,

S.B.

Jiang ,

C.M.

Ma ,

G.H.

Olivera ,

R.F.

Rodebaugh and

K.J.

Ruchala , The management of imaging dose during image-guided radiotherapy: Report of the AAPM Task Group 75, Medical Physics 34 (2007), 4041–4063.

D.J.

Brenner and

E.J.

Hall , Computed tomography—an increasing source of radiation exposure, New England Journal of Medicine 357 (2007), 2277–2284.

Feldkamp ,

Davis and

Kress , Practical cone-beam algorithm, JOSA A 1 (1984), 612–619.

Tang ,

Hsieh ,

R.A.

Nilsen ,

Dutta ,

Samsonov and

Hagiwara , A three-dimensional-weighted cone beam filtered backprojection (CB-FBP) algorithm for image reconstruction in volumetric CT-helical scanning, Physics in Medicine & Biology 51 (2006), 855.

Mori ,

Endo ,

Komatsu ,

Kandatsu ,

Yashiro , and

Baba , A combination-weighted Feldkamp-based reconstruction algorithm for cone-beam CT, Physics in Medicine & Biology 51 (2006), 3953.

10.

Wang ,

Li and

Xing , Iterative image reconstruction for CBCT using edge-preserving prior, Medical Physics 36 (2009), 252–260.

11.

E.Y.

Sidky and

Pan , Image reconstruction in circular cone-beam computed tomography by constrained, totalvariation minimization, Physics in Medicine & Biology 53 (2008), 4777–4807.

12.

Jia ,

Lou ,

Li ,

W.Y.

Song and

S.B.

Jiang , GPU-based fast cone beam CT reconstruction from undersampled and noisy projection data via total variation, Medical Physics 37 (2010), 1757–1760.

13.

Chen ,

and

Leng , Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Medical Physics 35 (2008), 660–663.

14.

Zheng ,

Cai ,

Li ,

Yan ,

Le ,

Wang and

Hu , Sparse-view image reconstruction via total absolute curvature combining total variation for x-ray computed tomography, Journal of X-ray Science and Technology 25 (2017), 959–980.

15.

Gao , Fast parallel algorithms for the x-ray transform and its adjoint, Medical Physics 39 (2012), 7110–7120.

16.

Han and

Yuan , A note on the alternating direction method of multipliers, Journal of Optimization Theory and Applications 155 (2012), 227–238.

17.

J.T.

Oden ,

Glowinski and

P.L.

Tallec , Augmented Lagrangian and Operator Splitting Method in Non-linear Mechanics, Mathematics of Computation 58 (1989).

18.

Glowinski and

J.T.

Oden , Numerical methods for nonlinear variational problems, Journal of Applied Mechanics 52 (1985), 739.

19.

Wang ,

A.C.

Bovik ,

H.R.

Sheikh and

E.P.

Simoncelli , Image quality assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing 13 (2004), 600–612.

20.

The Cancer Imaging Archive, http://www.cancerimagingarchive.net. Accessed 20 Dec 2016.

21.

Suyun ,

Ping ,

Jinxiao and

Yayu , A soft-threshold filtering EM-TV algorithm for CT reconstruction, Optik-International Journal for Light and Electron Optics 126 (2015), 202–207.

22.

Sawatzky ,

Brune ,

Wubbeling ,

Kosters ,

Schafers and

Burger , Accurate EM-TV algorithm in PET with low SNR, in: Nuclear Science Symposium Conference Record, 2008. NSS’08. IEEE, IEEE, 2008, pp. 5133–5137.

23.

Mascolo-Fortin ,

Matenine ,

Archambault and

Despres , A fast 4D cone beam CT reconstruction method based on the OSC-TV algorithm, Journal of X-ray Science and Technology 26 (2018), 189–208.

24.

Xu ,

Yang ,

Tan ,

Sawatzky and

M.A.

Anastasio , Accelerated fast iterative shrinkage thresholding algorithms for sparsity-regularized cone-beam CT image reconstruction, Medical Physics 43 (2016), 1849.