Performance of the sinogram-based iterative reconstruction in sparse view X-ray computed tomography

Abstract

Performing X-ray computed tomography (CT) examinations with less radiation has recently received increasing interest: in medical imaging this means less (potentially harmful) radiation for the patient; in non-destructive testing of materials/objects such as testing jet engines, the reduction of the number of projection angles (which for large objects is in general high) leads to a substantial decreasing of the experiment time. In the experiment, less radiation is usually achieved by either (1) reducing the radiation dose used at each projection angle or (2) using sparse view X-ray CT, which means significantly less projection angles are used during the examination. In this work, we study the performance of the recently proposed sinogram-based iterative reconstruction algorithm in sparse view X-ray CT and show that it provides, in some cases, reconstruction accuracy better than that obtained by some of the Total Variation regularization techniques. The provided accuracy is obtained with computation times comparable to other techniques. An important feature of the sinogram-based iterative reconstruction algorithm is that it is simpler and without the many parameters specific to other techniques.

Keywords

Dose reduction iterative reconstruction sparse view X-ray CT

1 Introduction

Iterative reconstruction algorithms [1, 2] have been extensively applied recently for the special case of sparse view X-ray computed tomography (CT) [1]. Iterative algorithms have been used for a long time in tomography, but only during the last few years several manufacturers have made available and suggested the use of iterative techniques for medical CT imaging that simultaneously provide for good image quality with detectability of low-contrast lesions and significant dose reduction. In the special case of sparse view X-ray CT (that is, dose reduction by using significantly less projection angles, or views, during the examination) the iterative reconstruction algorithms are applied for obtaining a reconstruction of accuracy comparable with the case that uses an usual number of views.

Besides the optimization of the Algebraic Reconstruction Technique (ART), Simultaneous Algebraic Reconstruction Technique (SART), and Simultaneous Iterative Reconstruction Technique (SIRT) methods [3 –5], the recent researches based on Total Variation (TV) regularization propose some of the most promising algorithms for sparse view X-ray CT [6 –14]. However, all these papers report long computation times associated with TV regularization techniques. For example:

In [6], computation times are reported for reconstructions of size 256 by 256 pixels and 90 views; they are between 19 and 28 seconds.

In [12], it is stated regarding the computation time, that "Similar to most SIR methods in CT image reconstruction (Ma et al, 2012a, Lauzier ThÃ©riault and Chen, 2013), the computational cost of the present PWLS-TGV method is very large because of the projection and back-projection operations using a huge system matrix. However, with the development of fast computers and dedicated hardware (Xu and Mueller, 2005, Knaup et al, 2006), iterative reconstruction algorithms, including the present PWLS-TGV method, may be used in clinical CT reconstruction in the near future." The authors did not report any computation time.

In [13] the authors have proposed TV regularization techniques that are run with 100 iterations, and each iteration takes 44.33 seconds, for a reconstruction of size 512 by 512.

Besides this, most papers developing TV regularization techniques report problems with low contrast structures that tend to be smoothed out by the TV regularization, posing a great challenge for these techniques.

In [15], which is representative of another set of works based on POCS (projection onto convex sets) for sparse view X-ray CT, the authors report a number of comparisons but do not report any computation time; instead it is said: "computation time for any reconstruction algorithm and, in particular, for iterative algorithms can be of concern. Accelerated computation can be achieved by streamlining/parallelizing the algorithms and/or by exploiting the available, or rapidly available, high-performance commodity computational hardware such as multi-core CPU and graphic processing units".

Other important works include: the works of Fessler et al. [16 –18] that propose some methods for accelerating the iterative reconstructions with focus on statistical iterative reconstruction; the works of Bouman et al. [19, 20] that propose fast methods for model-driven X-ray CT; and the works of Unser et. al that propose second-order extensions of TV regularization [21].

In [22] the performance of the recently proposed sinogram-based iterative reconstruction algorithm [23] was briefly studied on a Shepp-Logan phantom. In this work we study in more detail the performance of proposed sinogram-based iterative reconstruction algorithm, with focus on the case of sparse view X-ray CT.

It is shown that the sinogram-based iterative reconstruction algorithm provides good reconstruction accuracy (in some cases obtaining accuracy slightly better than that obtained by some variants of the TV regularization technique), with computation times comparable with those reported for the TV regularization techniques. The examples studied are: the NCAT thoracic cross-section studied also in other recent works, and then an image showing a pipe welding process cross-section. The important feature of the proposed algorithm is its simplicity in comparison with other techniques, in particular being much simpler and without the many parameters specific to the TV regularization techniques.

2 Sinogram-based iterative reconstruction

2.1 Review

First, we review the sinogram-based iterative reconstruction (SbIR) algorithm [22]. Let μ₁ be the density matrix (of size n_x lines by n_y columns) to be reconstructed from sinogram y, n_d the number of detectors and n_v the number of views. If M is n_d * n_v and N is n_x * n_y, then the system matrix A is of size M lines by N columns and the sinogram y is a vector of size M. If we write μ₁ in vector format then μ₁ is a vector of size N; the reconstruction problem is to find μ₁ such that

$A * μ_{1} = y .$

There are two steps in the reconstruction process of the SbIR algorithm:

initialization: first we initialize the reconstruction μ₁ with an approximation derived from the sinogram y as follows:

μ₁ ← A′ * (y./ α)

μ₁ ← μ₁./ β

iterations: at each iteration, we update the reconstruction μ₁ based on the difference between the sinogram y and the current situation in μ₁:

$\tilde{y}$ ← A * μ₁

μ₂ ← diag(μ₁) * A′ * (y./ $\tilde{y}$ )

μ₁ ← μ₂./ β

where:

A′ is the transpose of A,

./ is the element-by-element division,

$\tilde{y}$ is the situation in the reconstruction μ₁ at the current iteration, and is a vector of size M as the sinogram y,

α = ∑_jA [i, j] is of size M,

β = ∑_iA [i, j] is of size N;

each element of μ₁ corresponds to a square of size 1.0 by 1.0 in the experiment;

we also consider that y has been obtained as area integrals, at a Detector, View pair meaning that each of the elements that are traversed is traversed for an area ≤ 1.0; we could consider line integrals as well, but using area integrals seem to help in the understanding of the algorithm.

The algorithm is very natural and its workings can be explained as follows; first let’s explain the initialization (when looking at the steps 1.1, 1.2 above, look only at the first element of μ₁ which is μ₁ [1]):

in step 1.1, each of the M elements of y./ α is the average sinogram measurement per unit of area traversed by the X-ray beam at the respective Detector, View pair; thus each element of y./ α can be viewed as an approximation of each of the elements of μ₁ that are traversed by the X-ray beam at the respective Detector, View pair;

so, in step 1.1 the first element of A′ * (y./ α) is the sum of all approximations of μ₁ [1] (from all Detector, View pairs in which μ₁ [1] is traversed);

in step 1.2 μ₁ [1] is divided by β [1], where β [1] is the sum of all non-zero elements in the first line of A′; this is natural, because each of the approximations of μ₁ [1] that come from step 1.1 are for unit of area 1.0, but for some Detector, View pairs in which μ₁ [1] is traversed it might be traversed only for some portion of its total area 1.0.

Given this explanation of the initialization step, a better writing of the algorithm would be:

initialization: first we initialize the reconstruction μ₁ with an approximation derived from the sinogram y as follows:

μ₁ ← B * (y./ α)

iterations: at each iteration, we update the reconstruction μ₁ based on the difference between the sinogram y and the current situation in μ₁:

$\tilde{y}$ ← A * μ₁

μ₁ ← diag(μ₁) * B * (y./ $\tilde{y}$ )

where:

B is a matrix of size N lines by M columns, whose i-th line is obtained from the i-th line of A′ by dividing each element of the i-th line of A′ to β [i].

Given this re-writing of the algorithm, the workings can be explained as follows (again look only at μ₁ [1]):

in step 1.1, μ₁ [1] is obtained (as a weighted sum) by multiplying the first line of B with the vector of approximations y./ α; precisely, consider that there are k (k ≤ M) non-zero elements in the first line of B, let these be B [1, j₁], B [1, j₂],…, B [1, j_k] where 1 ≤ j₁ ≤ j₂ ≤ … ≤ j_k ≤ M; each of B [1, j₁], B [1, j₂],…, B [1, j_k] corresponds to one of those Detector, View pairs that traverse μ₁ [1], and represents the fraction (weight) from the total sum traversed β [1]; so we have

$B [1, j_{1}] + B [1, j_{2}] + . . . + B [1, j_{k}] = 1.0;$

thus if in step 1.1 above Y_α is the vector of approximations y./ α, then we have

$μ_{1} [1] \leftarrow B [1, j_{1}] * Y_{α} [j_{1}] + B [1, j_{2}] * Y_{α} [j_{2}] + . . . + B [1, j_{k}] * Y_{α} [j_{k}],$

which means that μ₁ [1] is obtained as a weighted sum of approximations;

the same explanation goes for the steps 2.1 and 2.2 from each iteration: in 2.1, $\tilde{y}$ is the situation of the current iteration, and thus the vector y./ $\tilde{y}$ is the vector of correction coefficients;

to see how μ₁ [1] is computed in 2.2, observe that it is computed by multiplying the first line of the matrix diag(μ₁) * B with the vector of corrections y./ $\tilde{y}$ ; thus if $Y_{\tilde{y}}$ is the vector of corrections y./ $\tilde{y}$ , then

$μ_{1} [1] \leftarrow B [1, j_{1}] * μ_{1} [1] * Y_{\tilde{y}} [j_{1}] + B [1, j_{2}] * μ_{1} [1] * Y_{\tilde{y}} [j_{2}] + . . . + B [1, j_{k}] * μ_{1} [1] * Y_{\tilde{y}} [j_{k}],$ which means that at each iteration μ₁ [1] is obtained as a weighted sum of approximations (the weights being again B [1, j₁], B [1, j₂],…, B [1, j_k]).

2.2 A simple example

Let n_x = n_y = n_d = n_v = 2 (so N = M = 4) and $μ_{1} = [\begin{matrix} 1.0 & 3.0 \\ 2.0 & 4.0 \end{matrix}]$ the matrix to be reconstructed, written in vector format as $μ_{1} = [\begin{matrix} 1.0 \\ 2.0 \\ 3.0 \\ 4.0 \end{matrix}] .$

We consider that the sinogram y is obtained without any noise using the views V1 and V2, the view V1 at 0 degrees and the view V2 at a rotation around the object of 90 degrees, as in Figure 1; if each of the 4 squares of μ₁ is of size 1.0 by 1.0, and the distance from X-ray source S to center of rotation (i.e. the center of reconstruction μ₁) is 2.0, and the distance from center of rotation to the detectors’ line is also 2.0, and the detectors’ line is of length 4.0, this means that:

in view V1, μ₁ [1] and μ₁ [2] are each traversed by a full area of 1.0, whereas μ₁ [3] and μ₁ [4] are each traversed by an area of 0.75;

in view V2, μ₁ [2] and μ₁ [4] are each traversed by a full area of 1.0, whereas μ₁ [1] and μ₁ [3] are each traversed by an area of 0.75.

Fig. 1

At the left is the geometry of view V1, and at the right the geometry of view V2; in yellow is the area of μ₁ that is traversed by that portion of the beam corresponding to detector D1 and in green is the area of μ₁ that is traversed by that portion of the beam corresponding to detector D2.

Thus

$A = [\begin{matrix} 1.0 & 0.0 & 0.75 & 0.0 \\ 0.0 & 1.0 & 0.0 & 0.75 \\ 0.75 & 1.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.75 & 1.0 \end{matrix}], y = [\begin{matrix} 3.25 \\ 5.0 \\ 2.75 \\ 6.25 \end{matrix}], α = [\begin{matrix} 1.75 \\ 1.75 \\ 1.75 \\ 1.75 \end{matrix}], β = [\begin{matrix} 1.75 \\ 2.0 \\ 1.5 \\ 1.75 \end{matrix}], B = [\begin{matrix} \frac{1.0}{1.75} & 0.0 & \frac{0.75}{1.75} & 0.0 \\ 0.0 & \frac{1.0}{2.0} & \frac{1.0}{2.0} & 0.0 \\ \frac{0.75}{1.5} & 0.0 & 0.0 & \frac{0.75}{1.5} \\ 0.0 & \frac{0.75}{1.75} & 0.0 & \frac{1.0}{1.75} \end{matrix}];$ the first line of A corresponds to D1,V1, the second line to D2,V1, the third line to D1,V2, and the fourth line to D2,V2; also note that the fact that all values in α are equal is particular of this example (in other geometries they could be all different).

The SbIR algorithm works as follows (look only at μ₁ [1]):

in 1.1, μ₁ [1] is initialized by multiplying the first line of B with the vector y./ α, so $μ_{1} [1] = \frac{1.0}{1.75} * \frac{3.25}{1.75} + \frac{0.75}{1.75} * \frac{2.75}{1.75} = 1.734;$ note that $\frac{3.25}{1.75}$ is the average sinogram measurement per unit of area in the D1,V1 pair and $\frac{2.75}{1.75}$ is the average sinogram measurement per unit of area in the D1,V2 pair; also note that μ₁ [1] is traversed in all Detector, View pairs for a total area of 1.75 (this is β [1]), with an area of 1.0 traversed in the D1,V1 pair and an area of 0.75 in the D1,V2 pair;

in 2.1, in the first iteration, the situation is $\tilde{y} = [\begin{matrix} 3.769 \\ 4.662 \\ 3.514 \\ 5.3 \end{matrix}] .$

in 2.2, in the first iteration, μ₁ [1] is obtained by multiplying the first line of diag(μ₁) * B with the vector $Y_{\tilde{y}}$ : $μ_{1} [1] = \frac{1.0}{1.75} * 1.734 * \frac{3.25}{3.769} + \frac{0.75}{1.75} * 1.734 * \frac{2.75}{3.514} = 1.434;$ before this first iteration the situation in μ₁ corresponding to the D1,V1 pair is 3.769, this means the elements that are traversed in D1,V1 need to be corrected by multiplication with $\frac{3.25}{3.769}$ , so this is where the first approximation $1.734 * \frac{3.25}{3.769}$ of μ₁ [1] comes from; before this first iteration the situation in μ₁ corresponding to the D1,V2 pair is 3.514, this means the elements that are traversed in D1,V2 need to be corrected by multiplication with $\frac{2.75}{3.514}$ , so this is where the second approximation $1.734 * \frac{2.75}{3.514}$ of μ₁ [1] comes from;

as we can observe, μ₁ [1] converges to its real value 1.0;

3 Performance of the SbIR algorithm, with focus on sparse view X-ray CT

We examine the performance of the SbIR algorithm by two examples, using the following scanning geometry for a fan-beam X-ray source, with detectors aligned in line:

the number of detectors is n_d = 1024, the size of each detector being 1.41, so the detectors line has length 1448.15;

the distance from the X-ray source to the center of rotation is 1024.0, and the distance from the center of rotation to the detectors’ line is also 1024.0;

the reconstruction is of size 512 by 512 pixels, with each pixel corresponding to a physical square area of size 1.0 by 1.0;

the number of views is n_v = 120 covering the full circle, so the angle (in radians) between two consecutive views is $\frac{2 π}{120}$ .

The SbIR algorithm is compared against the TV regularization techniques GPBB, UPN and UPN0 implemented in the TVReg software package developed recently at the Technical University of Denmark (DTU):

the Barzilai-Borwein accelerated gradient projection first-order method (GPBB);

implementation of an optimal first-order method for strongly convex functions, due to Nesterov, tailored to large-scale total variation regularization (UPN);

a first-order method by Nesterov, Beck, and Teboulle (UPN0) for the case of a zero strong convexity parameter.

The TVReg software package is written in C language, and is available at the Github software repository [25] and described in [26]. The TVReg package solves the problem $\min_{x} {\frac{1}{2} | | Ax - y | |_{2}^{2} + α TV (x)},$ where x represents the reconstructed image, y the sinogram data, and A the forward operator (the system matrix). TV (x) is the Total Variation of the image, and α is a regularization parameter set by the user. Other parameter (besides α) that can be set by the user is the noise parameter rnl that is used to add Gaussian noise to the synthetically constructed sinogram y in the case of simulated experiments.

Regarding the SbIR algorithm, the synthetically simulated sinograms in the two examples that we show are calculated using line integrals, as in most works in sparse view X-ray CT.

3.1 Example 1

In the first example we consider the NCAT phantom studied also in other works [7] and shown in Fig. 2(a) (this phantom simulates a thoracic cross-section); in this example we add noise to the sinogram by considering rnl = 0.005 for both SbIR and the GPBB, UPN, UPN0 algorithms. Also the number of iterations is 512 for all algorithms (SbIR, GPBB, UPN, UPN0).

The reconstruction with the SbIR algorithm is shown in Fig. 2(c), where in Fig. 2(b) is shown the reconstruction obtained after the initialization step is applied. The reconstructions with the GPBB, UPN, UPN0 algorithms are shown in Fig. 3 when the regularization parameter α is 0.5, 5.0, or 50.0. In this example, the best reconstructions obtained with the GPBB, UPN, UPN0 algorithms are obtained when α = 5.0 (this means the middle row (b), (e) and (h) in Fig. 3

Fig. 2

(a) The original NCAT phantom, (b) the initialization of the SbIR algorithm and (c) the reconstruction with the SbIR algorithm after 512 iterations; both (b) and (c) are obtained after adding noise (with rnl = 0.005) to the sinogram.

Fig. 3

(a), (b) and (c) are the results after 512 iterations for the GPBB method for α = 0.5, α = 5.0 and α = 50.0 respectively; (d), (e) and (f) are the results after 512 iterations for the UPN method for α = 0.5, α = 5.0 and α = 50.0 respectively; (g), (h) and (i) are the results after 512 iterations for the UPN0 method for α = 0.5, α = 5.0 and α = 50.0 respectively.

We calculated the SSIM (Structural Similarity Index), PSNR (Peak Signal-to-Noise Ratio), MSE (Mean Square Error) for all the reconstructions: high values of SSIM and PSNR means higher accuracy while high values of MSE means lower accuracy; for the SbIR reconstruction they are SSIM:0.9018, PSNR: 25.6149, MSE: 0.0027; for GPBB, UPN, UPN0 they are listed in Table 1; as we can see the only case when the reconstruction with the SbIR algorithm is better is the one shown in Fig. 3(c).

Table 1

calculations showing the difference between each of the reconstructions with GPBB, UPN, UPN0 and the original, for example 1

	GPBB	UPN	UPN0
α = 0.5	SSIM: 0.9427, PSNR: 37.1377, MSE: 1.9330e-04	SSIM: 0.9444, PSNR: 38.4690, MSE: 1.4226e-04	SSIM: 0.9447, PSNR: 38.4914, MSE: 1.4154e-04
α = 5.0	SSIM: 0.9791, PSNR: 35.1269, MSE: 3.0712e-04	SSIM: 0.9920, PSNR: 41.8835, MSE: 6.4812e-05	SSIM: 0.9919, PSNR: 41.9283, MSE: 6.4146e-05
α = 50.0	SSIM: 0.8869, PSNR: 27.3581, MSE: 0.0018	SSIM: 0.9690, PSNR: 32.1766, MSE: 6.0582e-04	SSIM: 0.9725, PSNR: 32.6921, MSE: 5.3801e-04

3.2 Example 2

A problem with the GPBB, UPN, UPN0 algorithms is that for a low value of α the reconstruction is of low quality, whereas for a high value of α the reconstruction tends to become "foggy" and with visible patches; in example 1, the value α = 5.0 seems to be a good choice.

We applied the same experimental setup but this time with rnl = 0.001 (so with lower level of noise) to the abdominal phantom shown in Fig. 4, proposed by the Institute of Medical Physics, Friedrich-Alexander-University Erlangen-Nurnberg, Erlangen [24]; the reconstruction with UPN for this phantom was also best for α = 5.0 as in the case of example 1, whereas GPBB and UPN0 were best for α = 0.5 in this case. The conclusion would be that α = 5.0 works best in general for the UPN algorithm, for reconstructions of size 512 by 512 pixels; for GPBB and UPN0, the choice of α is more unclear.

Fig. 4

Abdominal phantom proposed by the Institute of Medical Physics, Friedrich-Alexander-University Erlangen-Nurnberg, Erlangen.

In this second example we show that setting α = 5.0 could lead to reconstruction obtained by the UPN algorithm of lower quality than the one obtained with the SbIR algorithm. Consider the image shown in Fig. 5 (a), which corresponds to a cross-section from a pipe welding process, image acquired during 2015 at the Institute of Non-destructive Testing, Tomsk Polytechnic University, Tomsk (this is an image from a non-destructive material testing experiment instead of medical imaging, but it is equally important since reducing the number of views yields an overall reduction of the experiment time); in this example, we add noise to the sinogram by considering rnl = 0.001 for both SbIR and the GPBB, UPN, UPN0 algorithms. Also the number of iterations is 512 for all algorithms (SbIR, GPBB, UPN, UPN0).

The reconstruction with the SbIR algorithm is shown in Fig. 5(c), where in Fig. 5(b) is shown the reconstruction obtained after the initialization step is applied. The reconstructions with the GPBB, UPN, UPN0 algorithms are shown in Fig. 6 when the regularization parameter α is 0.5, 5.0, or 50.0. The best reconstructions obtained with the GPBB, UPN, UPN0 algorithms are obtained when α = 0.5 (this means the first row (a), (d) and (g) in Fig. 6).

Fig. 5

(a) The original pipe welding process image, (b) the initialization of the SbIR algorithm and (c) the reconstruction with the SbIR algorithm after 512 iterations; both (b) and (c) are obtained after adding noise (with rnl = 0.001) to the sinogram.

Fig. 6

We calculated the SSIM, PSNR, MSE for all the reconstructions; for the SbIR reconstruction they are SSIM: 0.8820, PSNR: 35.5374, MSE: 2.7942e-04; for GPBB, UPN, UPN0 they are listed in Table 2 as we can see the reconstruction with the SbIR algorithm is better than all 6 reconstructions from Fig. 6(b), (c), (e), (f), (h), (i), but worse than the reconstructions corresponding to α = 0.5 shown in Fig. 6(a), (d), (g). Thus, this second example shows that even if in general the TV regularization techniques converge faster than the SbIR algorithm, there are cases as the one shown here where the SbIR algorithm provides a better reconstruction. Also, the SbIR algorithm has no parameters to be set as in the case of the TV regularization techniques.

Table 2

calculations showing the difference between each of the reconstructions with GPBB, UPN, UPN0 and the original, for example 2

	GPBB	UPN	UPN0
α = 0.5	SSIM: 0.9077, PSNR: 37.4524, MSE: 1.7979e-04	SSIM: 0.9075, PSNR: 37.5485, MSE: 1.7585e-04	SSIM: 0.9066, PSNR: 37.5471, MSE: 1.7591e-04
α = 5.0	SSIM: 0.8760, PSNR: 34.5804, MSE: 3.4831e-04	SSIM: 0.8790, PSNR: 35.1570, MSE: 3.0500e-04	SSIM: 0.8790, PSNR: 35.1281, MSE: 3.0703e-04
α = 50.0	SSIM: 0.7947, PSNR: 28.4132, MSE: 0.0014	SSIM: 0.8372, PSNR: 31.7535, MSE: 6.6780e-04	SSIM: 0.8382, PSNR: 31.8891, MSE: 6.4728e-04

Fig. 7

Convergence of the SbIR algorithm, example 2 (the pipe welding process image).

3.3 Computation times achieved by the SbIR, GPBB, UPN, UPN0 algorithms

The reconstructions have been obtained on a conventional PC with an i7-6500u @ 2.5 GHz 4-cores processor, and 8 GB of memory; only one core out of the four available has been used for each of the algorithms SbIR, GPBB, UPN, UPN0. The SbIR algorithm was implemented in C language in Visual C++ 2015; the GPBB, UPN, UPN0 algorithms of the TVReg package are implemented in C language. For the reconstructions shown in the two examples (all with 512 iterations) the computation time was around 150 seconds for the SbIR algorithm and between around 100 and 250 seconds for the GPBB, UPN, UPN0 algorithms.

3.4 Convergence of the SbIR algorithm

The convergence of the SbIR algorithm, after a number of studies we have done, seems to behave as follows: first the SbIR algorithm spends a number (that depends on the specific input and noise in the sinogram) of iterations to stabilize within a certain interval, after which the rest of the iterations are spent to converge to the optimal solution. For the second example (the pipe welding process image) that is run with 512 iterations, the convergence behavior is shown in Fig. 2 the red graph shows the sum $\sum_{i = 1}^{M} y [i]$ of all sinogram values, which for the second example it is 8.987062 × 10⁶; and in blue graph it is plotted the sum $\sum_{i = 1}^{M} \tilde{y}$ of all situation values, at each of the 512 iterations. We observe that in this plot the SbIR algorithm spends around 100 iterations to stabilize within a certain interval, and after that starts to converge to an approximate solution.

We have tried other initial solutions for the SbIR algorithm (like the solution given by the Filtered Back-Projection algorithm) to see if with other initializations it converges faster, but it seems that its initialization as shown in this work works best.

4 Conclusions

From the results reported, the SbIR algorithm clearly proves to be a good option for obtaining good reconstruction accuracy in the special case of sparse view X-ray CT, in some cases with accuracy slightly better than some variants of the TV regularization technique. Comparing to the TV regularization techniques the SbIR algorithm is much simpler and without any parameter. One important work that needs to be done is proving the convergence of the SbIR algorithm mathematically. Other important further works include: (1) trying to find better initial solutions in order to speed up its convergence; and (2) development of optimal parallelizations using either multi-core processors, or Graphical Processing Units (GPUs), or both (hybrid CPU-GPU parallelizations).

References

Beister ,

Kolditz and

W.A.

Kalender , Iterative reconstruction methods in X-ray CT, Physica Medica: European Journal of Medical Physics 28 (2012), 94–108.

A.C.

Kak and

Slaney , Principles of computerized tomographic imaging, Series: Classics in Applied Mathematics (Book 33), Society for Industrial and Applied Mathematics, 2001.

Gordon ,

Bender and

G.T.

Herman , Algebraic reconstruction techniques (ART) for threedimensional electron microscopy and x-ray photography, Journal of Theoretical Biology 29 (1970), 471–482.

A.H.

Andersen and

A.C.

Kak , Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm, Ultrasonic Imaging 6 (1984), 81–94.

Gilbert , Iterative methods for the three-dimensional reconstruction of an object from projections, Journal of Theoretical Biology 36 (1972), 105–117.

Zhang ,

Wang ,

Yan ,

Li ,

Cai and

Hu , Constrained total generalized p-variation minimization for few-view X-ray computed tomography image reconstruction, PLoS One 11 (2016), paper e0149899, pp. 28.

Tian ,

Jia ,

K.H.

Yuan ,

Pan and

S.B.

Jiang , Low-dose CT reconstruction via edge-preserving total variation regularization, Physics in Medicine and Biology 56 (2011), 5949–5967.

Liu ,

J.H.

Ma ,

Fan and

Z.R.

Liang , Adaptive-weighted total variation minimization for sparse data toward low-dose X-ray computed tomography image reconstruction, Physics in Medicine and Biology 57 (2012), 7923–7956.

Cai ,

Wang ,

Zhang ,

Yan ,

Li ,

Xi and

Li , Edge guided image reconstruction in linear scan CT by weighted alternating direction TV minimization, Journal of X-Ray Science and Technology 22 (2014), 335–349.

10.

Chang ,

Li ,

Chen ,

Xiao ,

Zhang and

Wang , A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction, Journal of X-Ray Science and Technology 21 (2013), 161–176.

11.

J.S.

Yang ,

H.Y.

Yu ,

Jiang and

Wang , High-order total variation minimization for interior tomography, Inverse Problems 26 (2010), 035013.

12.

Niu ,

Gao ,

Bian ,

Huang ,

Chen ,

Yu ,

Liang and

Ma , Sparse-view X-ray CT reconstruction via total generalized variation regularization, Physics in Medicine and Biology 59 (2014), 2997–3017.

13.

Qi ,

Chen ,

Guo and

Zhou , Sparse-view computed tomography image reconstruction via a combination of L1 and SL0 regularization, Bio-Medical Materials and Engineering 26 (2015), S1389–S1398.

14.

Li ,

Zhang ,

Peng ,

Guan ,

Xu ,

Sun and

Zheng , Smoothed l0 norm regularization for sparse-view X-ray CT reconstruction, Biomedical Research International 2016 (2016), paper 2180457, 12.

15.

Han ,

Bian ,

D.R.

Eaker ,

T.L.

Kline ,

E.Y.

Sidky ,

E.L.

Ritman and

Pan : Algorithm-enabled low-dose micro-CT imaging, IEEE Transactions on Medical Imaging 30 (2011), 606–620.

16.

Ramani and

J.A.

Fessler , A splitting-based iterative algorithm for accelerated statistical X-ray CT reconstruction, IEEE Transactions on Medical Imaging 31 (2012), 677–688.

17.

J.K.

Kim ,

J.A.

Fessler and

Zhang , Forward-projection architecture for fast iterative image reconstruction in X-ray CT, IEEE Transactions on Signal Processing 60 (2012), 5508–5518.

18.

Kim ,

Pal and

J.-B.

Thibault , Accelerating ordered subsets image reconstruction for X-ray CT using spatially nonuniform optimization transfer, IEEE Transactions on Medical Imaging 32 (2013), 1965–1978.

19.

Yu ,

J.-B.

Thibault ,

C.A.

Bouman ,

K.D.

Sauer and

Hsieh , Fast model-based X-ray CT reconstruction using spatially non-homogeneous ICD optimization, IEEE Transactions on Image Processing 20 (2011), 161–175.

20.

Sabne ,

Wang ,

Kisner ,

Bouman ,

Raghunathan and

Midkiff , Model-based iterative CT image reconstruction on GPUs, 21st ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (PPoPP'17) (2017), 207–220.

21.

Dogan ,

Lefkimmiatis ,

Bourquard and

Unser , A second-order extension of TV regularization for image deblurring, Proceedings of the IEEE International Conference on Image Processing (2011), 713–716.

22.

Trinca ,

Zhong ,

Y.-Z.

Wang ,

Mamyrbayev and

Libin , Sinogram-based adaptive iterative reconstruction for sparse view X-ray computed tomography, Proceedings of SPIE/COS Photonics Asia, vol. 10020 (2016), paper 100201G, pp. 6.

23.

D. Trinca, Implementation in Visual C++ 2015 of the sinogram-based iterative reconstruction, available at the Github software repository. https://github.com/dntrinca

24.

Phantom database, Institute of Medical Physics, Friedrich-Alexander-University Erlangen-Nurnberg, Erlangen, Germany; http://www.imp.uni-erlangen.de/phantoms/index.htm

25.

TVReg: Software package for 3D Total Variation regularization, available at the Github software repository https://github.com/jakobsj/TVReg

26.

T.L.

Jensen ,

J.H.

Jørgensen ,

P.C.

Hansen and

S.H.

Jensen , Implementation of an optimal first-order method for strongly convex total variation regularization, BIT Numerical Mathematics 52 (2012), 329–356.