1

Introduction

Linear models for computed tomography (CT) play an important role for iterative image reconstruction. The most common approach to CT processing involves taking the negative logarithm of the projection data, so that the line integration model leads to a linear relation between the image and processed data. Accordingly, the CT image reconstruction problem can be written as a large linear system

where b, a vector of lengthm, represents the processed projection data; x, a vector of length n, contains the image pixel values; and the m × n system matrix A contain the weights that model line-integration. Linear tomographic models can include quadratic regularization, cf.¹ and [2, Chapter 12], or more sophisticated modeling such as noise correlation and blur due to accurate detector physics.³ Even when non-linear models for CT⁴ are considered for iterative image reconstruction, there is usually a large linear system that is involved in the algorithm. Novel techniques for solving large linear systems may thus be of practical use for iterative image reconstruction in CT.

The most common iterative algorithm for solving linear CT models, excluding row-action, sequential, or SIRT-type data processing methods, is the conjugate gradients (CG) algorithm [2, Chapter 11]. For least-squares problems with non-symmetric system matrices, in particular, there is the conjugate-gradients least-squares (CGLS) algorithm that solves the optimization problem

The minimizer of this optimization problem can also be found from solving the normal equations directly

which are derived from (2) by taking the gradient of the objective function and setting it to zero. In applying CGLS, the implementation for back-projection B must be the matrix transpose A^⊤. If B ≠ A^⊤ then the resulting method is not well defined, there is no convergence theory, and if it does converge it does not solve Eqs. (2) and (3). Never-theless there are practical motivations for considering back-projection implementations B different than A^⊤. These motivations are outlined in Ref.⁵ in connection with SIRT-type iterative solvers, where the authors explain that B can be a preconditioner, B may be an efficient but approximate implementation of A^⊤, or A may involve complex physics modeling that may make computer implementation of A^⊤prohibitively expensive. As shown in⁶ we can guarantee convergence of SIRT-type methods with B ≠ A^⊤ (but not CGLS) by shifting the complex eigenvalue spectrum of BA so that eigenvalues with negative real part are eliminated; but it forces a modification of the problem that is being solved.

Use of the GMRES algorithm allows for use of back-projectors B that are not equal to A^⊤ without modification of the desired reconstruction model. Furthermore, the algorithm does not involve any parameters other than the iteration number. In Sec. 2 we present the ABBA framework⁷ which involves two forms of GMRES called AB-GMRES and BA-GMRES. In Sec. 3 we demonstrate use of BA-GMRES for unmatched projector/back-projector pairs and for preconditioning. We conclude this abstract in Sec. 4.

2

The ABBA framework

The GMRES algorithm solves a linear system

where the coefficient matrix S is a square matrix that is not necessarily symmetric. The relevance for CT image reconstruction is that a square non-symmetric matrix arises when multiplying unmatched back-projection B and projection A matrices; i.e., both AB and BA are square non-symmetric matrices. The original linear system of interest, Eq. (1), cannot be directly solved with GMRES because A is not necessarily square for CT, but this equation can be modified to

where the unknown vector y has the same length as the projection data b and the resulting m×m matrix AB is square. Also, the normal equations in Eq. (3) can be modified by replacing A^⊤ with B

and again the resulting n × n matrix BA is square. See⁸ for details. Modeling CT with Eq. (4) is similar to the use of natural pixels^9–11 as the image is expressed as the back-projection of a sinogram.

We refer to the GMRES algorithms for solving Eqs. (4) and (5) as AB-GMRES and BA-GMRES, respectively. The pseudo-code for both algorithms is given in Ref.,⁷ and we briefly describe the algorithms here. Similar to CGLS, GMRES is a Krylov subspace method, where the basis vectors of the subspace are generated by choosing an initial vector and repeatedly applying the coefficient matrix (AB or BA) to obtain new linearly independent basis vectors. For AB-GMRES or BA-GMRES with a zero initial vector, the first basis vector is b or Bb, respectively, and subsequent basis vectors are generated by applying the matrix AB or BA, respectively. The GMRES algorithm involves orthonormalization of the Krylov subspace vectors to obtain a orthonormal basis set that spans the subspace. At each GMRES iteration the dimension of the subspace is increased by one and the minimum residual that can be expressed by the Krylov basis set is found.

The computational burden of GMRES lies with the fact that the Krylov basis set must be stored and the number of basis vectors is the same as the number iterations. For AB-GMRES and BA-GMRES the size of one basis vector is the same as the size of a sinogram and image, respectively. In our implementation, the basis set is stored on the computer disk. Restart methods¹² can reduce the basis vector storage burden, but for this work we demonstrate the basic GMRES implementation.

The AB-GMRES and BA-GMRES algorithms are guaranteed to minimize different data discrepancy measures. In the case of AB-GMRES, the algorithm minimizes

while BA-GMRES minimizes

Note that BA-GMRES is not necessarily minimizing

In this work we focus on BA-GMRES and we demonstrate its use on cone-beam CT image reconstruction.

3

BA-GMRES applied to cone-beam CT image reconstruction

We apply BA-GMRES to a cone-beam CT (CBCT) data set acquired on an Epica Pegaso veterinary CT scanner. The particular scan configuration for the data set is 180 projections taken uniformly over one circular rotation. The detector size is 1088×896 detector pixels, where each pixel is (0.278mm)² in size. The 180-view dataset is sub-sampled from a 720-view scan of a quality assurance (QA) phantom. Image volumes are reconstructed onto a 1024 × 1024 × 300 voxel grid, and a reference volume is generated by use of filtered back-projection (FBP) applied to the full 720-view dataset and shown in Fig. 1. Also shown in the figure is FBP applied to the 180-view sub-sampled dataset.

Fig 1

Mid-slice images of QA phantom. (Left) FBP reconstructed image from 720 views. (Right) FBP reconstructed images from 180 views, a 4-fold sub-sampling of the original CBCT dataset. The grayscale window is [0.,0.25] cm⁻¹.

To demonstrate application of BA-GMRES to CBCT image reconstruction, we use a ray-driven cone-beam projector where the matrix elements for A are computed by the line-intersection method. We consider two implementations of B: (1) B_un voxel-driven back-projection using linear interpolation to determine the appropriate projection value on the detector, and (2) B_FBP = B_unF filtered back-projection, where F represents the ramp filter. See [2, Chapter 9] for details about these discretization models. The first BA-GMRES implementation tests unmatched back-projector/projector pairs where B_un ≈ A^⊤, and the second implementation includes the additional ramp-filtering step for preconditioning. Use of B_un and A is also shown for CGLS, which requires that B = A^⊤.

The data root-mean-square-error (RMSE) curves for both forms of BA-GMRES and CGLS using B_un and A are shown in the top panel of Fig. 2. The CGLS result initially shows convergence of the data RMSE, but after 20 iterations the data RMSE begins to diverge with increasing iteration number, as expected, since this algorithm is not designed to work for unmatched matrix transpose implementations. The corresponding BA-GMRES result does show a decreasing data RMSE with iteration number. For the preconditioned form of BA-GMRES, the decrease in data RMSE is even more rapid. The decreasing trends in ||Ax − b||₂ for BA-GMRES occur even though this algorithm is not guaranteed to reduce this data norm. Also shown in Fig. 2 is the data RMSE curves for ||BAx – Bb||₂, which is guaranteed to decrease with iteration number and they do indeed show decreasing trends for BA-GMRES. These issues are elaborated in.⁷

Fig 2

Data RMSE in the form of (Top) ||Ax − b||₂ and (Bottom) ||BAx − Bb||₂. In the top graph unmatched CGLS is also shown to demonstrate divergence of the data RMSE with unmatched projector/back-projector pairs. The other two curves correspond to use of voxel-driven back-projection В = B_un, and FBP В = B_FBP. The projector A is a ray-driven implementation.

The mid-slice images for BA-GMRES using both B implementations are shown in Fig. 3 at different iteration numbers. Preconditioning has a clear effect on the convergence as all the phantom structures are clearly visible in the early iterations and the gray-level is stabilized already at the fifth iteration. The BA-GMRES result without preconditioning is also fairly efficient as the main features of the QA phantom are visible at 20 iterations.

Fig 3

Mid-slice BA-GMRES images for voxel-driven back-projection B = B_un (Left column) and FBP B = B_FBP (Right column). The shown iteration numbers are 2, 5, 10, and 20 going from the top row to bottom row. The grayscale window is [0.,0.25] cm⁻¹.

One measure of image quality is to compare the reconstructed volumes to a ground truth image. Employing the 720-view FBP reconstructed volume as a surrogate for the ground truth, the image RMSE is plotted in Fig. 4 for both versions of BA-GMRES. The BA-GMRES implementation with B = B_un achieves a minimum image RMSE of 0.0201 at iteration 29, while the preconditioned version with B = B_un achieves a minimum image RMSE of 0.0210 at iteration 4. For comparison the 180-view FBP result has an image RMSE of 0.0347. To appreciate the various image qualities, ROI images of the mid-slice are shown at the minimum image RMSE iteration numbers in Fig. 5.

Fig 4

Using the 720-view reconstructed volume (see mid-slice image on the left of Fig. 1 as a reference), the BA-GMRES reconstructed image RMSE is plotted as a function of iteration number for B = B_un and B = B_FBP.

Fig 5

Mid-slice ROI images of QA phantom. (Top, Left) FBP reconstructed image from 720 views. (Top, Right) FBP reconstructed images from 180 views. (Bottom, Left) BA-GMRES image for B = B_un at iteration 29. (Bottom, Right) BA-GMRES image for B = B_FBP at iteration 4. The grayscale window is [0.,0.25] cm⁻¹.

That the image RMSE has a minimum at finite iteration number is a well-known phenomenon in iterative image reconstruction and it is known as semi-convergence [2, Chapter 11]. Early stopping in such algorithms is a form of regularization because the components associated with large singular values of A converge fast, while the unwanted noisy components associated with smaller singular values – that cause strong image artifacts – appear after more iterations. Semi-convergence is observed in the image RMSE curves of Fig. 4 and visually in the preconditioned BA-GMRES series of Fig. 3 where the image at 20 iterations clearly shows strong artifacts from iterating too far. The semi-convergence issue also presents a practical dilemma for preconditioning. With the shown preconditioned BA-GMRES results, the minimum image RMSE result is obtained already at the fourth iteration; thus the iteration number provides only coarse control over its image quality. The un-preconditioned BA-GMRES implementation achieves its image RMSE minimizer at the 29th iteration, which is computationally less efficient, but on the other hand the iteration number provides a finer control over the image quality. In any case, the BA-GMRES framework provides a flexible means for implementing back-projectors or preconditioning schemes, and optimizing the B implementation and iteration number will depend on the imaging task of interest.

4

Conclusion

This work presents an iterative image reconstruction framework for linear CT problems that allows for the use of unmatched back-projector/projector pairs in a straight-forwardmanner. This possibility is convenient for implementation of efficient back-projectors, linear modeling of complex physics, and preconditioning. Also, because it is clear what equation is being solved when B ≠ A^⊤, BA-GMRES can be used for solving linear sub-problems that may arise in non-linear iterative image reconstruction. The BA-GMRES algorithm does present a challenge for computer memory because the Krylov basis set needs to be stored during the iteration, but the present demonstration on CBCT image reconstruction does show that BA-GMRES can be applied to large-scale CT systems of clinical interest.