Singular value decomposition-based 2D image reconstruction for computed tomography

Abstract

Singular value decomposition (SVD)-based 2D image reconstruction methods are developed and evaluated for a broad class of inverse problems for which there are no analytical solutions. The proposed methods are fast and accurate for reconstructing images in a non-iterative fashion. The multi-resolution strategy is adopted to reduce the size of the system matrix to reconstruct large images using limited memory capacity. A modified high-contrast Shepp-Logan phantom, a low-contrast FORBILD head phantom, and a physical phantom are employed to evaluate the proposed methods with different system configurations. The results show that the SVD methods can accurately reconstruct images from standard scan and interior scan projections and that they outperform other benchmark methods. The general SVD method outperforms the other SVD methods. The truncated SVD and Tikhonov regularized SVD methods accurately reconstruct a region-of-interest (ROI) from an internal scan with a known sub-region inside the ROI. Furthermore, the SVD methods are much faster and more flexible than the benchmark algorithms, especially in the ROI reconstructions in our experiments.

Keywords

Computed tomography singular value decomposition a region of interest interior tomography Tikhonov regularization

1 Introduction

Although computed tomography (CT) has been widely applied for medical diagnosis, its potential radiation hazard cannot be ignored. CT scanning accounted for approximately 75% of the ionizing radiation exposure in 2002 in the United States [1]. Many methods have been proposed to reduce the excessive radiation dose, including three representative methods. The first is to reduce the number of photons emitted from the x-ray source; however, the noisy projections that are acquired degrade the image quality. The second method is to reduce the number of views. However, the reconstructed images contain severe streak artifacts. The third method is to illuminate a small region of interest (ROI) inside the object [2], which is referred to as an interior scan. However, due to the nonlocal property of the inverse Radon transform, no theoretical exact local reconstruction can be achieved solely from projections that pass through the ROI [3]. Based on prior knowledge, several theoretically exact solutions have been developed, which are referred to as interior tomography [2 , 4–8]. In the early 1990s, Gel’fand and Graev established a relationship between the Hilbert transform of an image along a line and the corresponding backprojection of the differentiated projections [9, 10]. The interior problem was proven to be uniquely solvable if there is a known sub-region inside the ROI [6, 7], which is referred to as knowledge-based interior tomography.

To date, two major reconstruction frameworks have been developed for knowledge-based interior tomography. One is an iterative framework that is based on projection onto convex sets [6]. An l₁ norm minimization of a sparse transform is usually applied as the regularization term. Iterative reconstruction algorithms involve multiple time-consuming projection and backprojection procedures, which limits their general applications. The other reconstruction framework is to use the singular value decomposition (SVD) method to reconstruct an image on 1D PI-segments based on the Hilbert transform relationship [11, 12]. Because this approach allows an interior reconstruction without iteration, it is much faster than its iterative counterparts. However, a rebinning procedure is usually required to convert the reconstructed images to Cartesian coordinates. Thus, the spatial resolution is compromised.

The greater capacity of computer memory motivates us to develop SVD-based 2D image reconstruction methods. These methods can perform fast image reconstructions for a broad class of inverse problems that have no analytical solutions. This includes, but is not limited to, few-view global reconstruction, local reconstruction from truncated projections, and irregular projection geometry. In particular, this paper emphasizes interior tomography as a practical application.

Several SVD-based methods have been proposed for single photon emission computed tomography (SPECT) reconstruction [13 –17]. Smith et al. [14] utilized a generalized matrix inverse to estimate the source activity distributions from SPECT. The SVD of the system matrix provided considerable insight into the SPECT image reconstruction problem. Tradeoffs between the resolution and error in estimating the source voxel intensities were also discussed in [14]. Zeng and Gullberg [15] applied SVD to analyze and reconstruct images. They concluded that problems might be ill-conditioned with truncated projections and that the reconstruction artifacts are due to the mismatch between the measured data and modeled projections. Wiener filtering and truncated inverse filtering were discussed in the search for optimal restoration with an SVD pseudo inversion method in [16]. Verhoeven [17] compared five CT algorithms for limited data problems, including SVD methods. SVD-based image reconstruction techniques have also been applied in other fields [18, 19]. However, a systematic evaluation of SVD methods for 2D CT image reconstruction has not been reported.

This paper is organized as follows. Section II presents mathematical models of CT imaging and SVD methods. The details of the algorithm implementation are presented in Section III, and Section IV shows the results. Finally, relevant issues are discussed in Section V.

2 Mathematical theory

2.1 Imaging model and SVD method

In a discrete CT system, both the image and the projection are modeled as vectors. The projection procedure is modeled as $Af = p,$ (1) where f ∈ ℝ^N is the image, p ∈ ℝ^M is the projection, and the projection procedure is simplified to a linear transform that is represented by a system matrix $A \in ℝ_{+}^{M \times N}$ , where ℝ is the real number domain, and ℝ ₊ represents a set of all nonnegative real numbers. An intuitive way to find the solution to Equation (1) is $\hat{f} {= A}^{p},$ (2) where A^† is the (pseudo) inverse of A. According to SVD, Acan be decomposed as $A = U Σ V^{T} \in ℝ_{+}^{M \times N},$ (3) where T represents the transpose of a matrix, bothU ∈ ℝ^M×M and V ∈ ℝ^N×N are unitary matrices, and ∑ ∈ ℝ^M×N is a uniquely determined diagonal matrix, whose diagonal entries σ_i (i ∈ { 1, 2, …, min(M, N) }) are singular values of A. The (pseudo) inverse of Acan be calculated as $A^{†} = V Σ^{†} U^{T},$ (4) where ∑^† is the (pseudo) inverse of ∑, which can be formed by replacing each non-zero diagonal entry σ_i by its reciprocal.

We can interpret the SVD of A in terms of eigen-decomposition. The columns of U are eigenvectors of AA^T, and the columns of V are eigenvectors of A^TA. In the CT reconstruction problem, if the column of U is backprojected and then projected, we can obtain the same projection at a constant scale. Using interior scan geometry to project a 128×128 image to 128 detector cells with 128 views, several representative columns of Uand Vare depicted as 2D images in Fig. 1. The first column ofU is the projection of a constant image, and the first column of Vis the backprojection of a constant projection data set.

2.2 Truncated SVD and Tikhonov regularization method for image reconstruction

Small singular values magnify projection noise. Therefore, severe artifacts will appear if the (pseudo) inverse is directly applied for image reconstructions. To suppress noise, we need to find a threshold to truncate small singular values. Given the magnitude ratio of singular values, the truncation position can be calculated according to $\hat{i} = \underset{i}{arg} min | \frac{\sum_{j = 1}^{i} σ_{j}}{\oplus σ} - r |,$ (5) where ⊕σ is the sum of all of the singular values, andr is a target ratio. If $ɛ (= σ_{\hat{i}})$ is a selected threshold, truncated SVD (TSVD) is usually applied to suppress the noise [20]. The pseudo-inverse of∑ is defined with ɛ as $σ_{i}^{- 1} = {\begin{matrix} 1 / σ_{i} & σ_{i} > ɛ \\ 0 & σ_{i} \leq ɛ \end{matrix} .$ (6)

A more general regularization method can be modeled as [20] ${∥ Af - p ∥}_{2}^{2} + ξ^{2} {∥ Φ f ∥}_{2}^{2},$ (7) where ξis a parameter that balances the data fidelity and regularization terms, and Φ represents a transform of f. If Φ = I, Equation (7) can be explicitly solved by defining the pseudo inverse of Σ $σ_{i}^{- 1} = \frac{σ_{i}}{σ_{i}^{2} + ξ^{2}} .$ (8)

Equation (8) is the solution of Tikhonov regularization [21], which is referred to as VSVD. Both TSVD and VSVD play similar roles in smoothing reconstructed images.

2.3 Generalized SVD method for image reconstruction

When Φ ≠ I, generalized SVD (GSVD [22]) can be used to solve Equation (7). Both A ∈ ℝ^M×N and Φ ∈ ℝ^P×N are real matrices with M ≥ N ≥ P and rank (Φ) = P. In GSVD, A and Φ can be decomposed as $\begin{array}{l} A = U (\begin{matrix} Λ & 0 \\ 0 & I \end{matrix}) B^{- 1} \\ Φ = V (\begin{matrix} M & 0 \\ 0 & 0 \end{matrix}) B^{- 1} \end{array}$ (9)

B is an invertible matrix. Both Λ and M are diagonal matrices, $\begin{matrix} Λ = diag (λ_{1}, λ_{2}, . . ., λ_{P}), 0 \leq λ_{1} \leq . . . \leq λ_{P} \leq 1, \\ M = diag (μ_{1}, μ_{2}, . . ., μ_{P}), 1 \geq μ_{1} \geq . . . \geq μ_{P} \geq 0 . \end{matrix}$ (10)

The values λ_i and μ_i are normalized so that $λ_{i}^{2} + μ_{i}^{2} = 1$ . Generalized singular values are defined as $γ_{i} = \frac{λ_{i}}{μ_{i}},$ (11) in a decreasing manner on the index i. A property of GSVD is $\begin{matrix} B^{T} A^{T} AB = (\begin{matrix} Λ^{2} & 0 \\ 0 & I \end{matrix}), \\ B^{T} Φ^{T} Φ B = (\begin{matrix} M^{2} & 0 \\ 0 & 0 \end{matrix}) . \end{matrix}$ (12)

The objective function Eq. (7) can be written in quadratic form as $\hat{f} = arg min_{x} {∥ (\begin{matrix} A \\ ξ Φ \end{matrix}) f - (\begin{matrix} p \\ 0 \end{matrix}) ∥}_{2}^{2}$ (13)

The optimal solution can be determined by setting the first derivative f with respect to zero, $(A^{T} A + ξ^{2} Φ^{T} Φ) f = A^{T} p .$ (14)

Therefore, $\hat{f} = {(A^{T} A + ξ^{2} Φ^{T} Φ)}^{- 1} A^{T} p .$ (15)

Equation (12) is applied to calculate (A^TA + ξ²Φ^TΦ) ^-1 with the following steps $\begin{matrix} B^{T} (A^{T} A + ξ^{2} Φ^{T} Φ) B = B^{T} A^{T} AB + ξ^{2} B^{T} Φ^{T} Φ B = (\begin{matrix} Λ^{2} + ξ^{2} M^{2} & 0 \\ 0 & I \end{matrix}), \\ (A^{T} A + ξ^{2} Φ^{T} Φ) = {(B^{T})}^{- 1} (\begin{matrix} Λ^{2} + ξ^{2} M^{2} & 0 \\ 0 & I \end{matrix}) B^{- 1}, \\ {(A^{T} A + ξ^{2} Φ^{T} Φ)}^{- 1} = {(B^{- 1})}^{T} {(\begin{matrix} Λ^{2} + ξ^{2} M^{2} & 0 \\ 0 & I \end{matrix})}^{- 1} B^{- 1} = {(B^{- 1})}^{T} \sum^{†} B^{- 1} . \end{matrix}$ (16) where Σ^† is a diagonal matrix that can easily be calculated. B is a full rank matrix, and its inverse needs to be calculated only once. In this paper, the Laplacian operator ∇ (= Φ) is applied.

2.4 Interior reconstruction with a known sub-region

Let Ω ={ 1, 2, …, N } be the index set of pixels. If we treat pixels as independent variables, A can be modified column by column. When a sub-region inside the ROI is known, the pixels can be grouped as $f = (\begin{matrix} f_{K} \\ f_{Ω ∖ K} \end{matrix}) \in ℝ^{N},$ (17) where f_K and f_Ω∖Krepresent unknown and known pixels, respectively, and K is the index set of known pixels. Correspondingly, the system matrix can be divided into two sub-matrices $A = [A_{K} A_{Ω ∖ K}] .$ (18)

The projection procedure is $Af = [A_{K} A_{Ω ∖ K}] (\begin{matrix} f_{K} \\ f_{Ω ∖ K} \end{matrix}) = A_{K} f_{K} + A_{Ω ∖ K} f_{Ω ∖ K} = p .$ (19)

We immediately arrive at the following linear matrix equation $A_{Ω ∖ K} f_{Ω ∖ K} = p - A_{K} f_{K} .$ (20)

Only the sub-matrix A_Ω∖K needs to be decomposed with SVD for the image reconstruction $A_{Ω ∖ K} = U_{Ω ∖ K} \sum_{Ω ∖ K} V_{Ω ∖ K}^{T}$ (21)

2.5 Multi-resolution image reconstruction with SVD methods

Reconstructing a high resolution image by SVD requires extensive memory resources. For example, to reconstruct a 512×512 image, A in a non-compressed format occupies 512 GB of double floating point data when the number of detector cells and the number of views are both 512. This challenging requirement inspires us to reconstruct the image using a multi-resolution scheme. In interior tomography, the ROI is strictly inside the object support. Therefore, using the multi-resolution strategy to compress the system matrix for ROI reconstruction will save space. We represent the object with fine pixels inside the ROI and with coarser pixels outside the ROI. Let A_c be an interior scan of a low resolution image and A_f be the same scan of a high resolution image. The columns of A_fthat correspond to the pixels inside the ROI are extracted to form A_L, and the columns of A_c that correspond to the pixels outside the ROI are extracted to formA_R. The interior scan system matrix of a dual-resolution image is constructed by merging A_L and A_R $A_{M} = [A_{L}, A_{R}] .$ (22)

The multi-resolution strategy provides a way to compress the system matrix. For example, if a 256×256 image is reconstructed using SVD methods, the size of A is 256² × 256² × 8 bytes = 32GB of double floating point data. Performing SVD requires at least 96 GB of memory. In Equation (22), the pixel size is 0.176 cm × 0.176 cm inside the ROI and 0.528 cm × 0.528 cm outside the ROI. Therefore, one coarse pixel covers nine fine pixels. A 153×153 fine grid is used to cover the ROI, and every 3×3 pixels outside the ROI is merged into one coarse pixel (see Fig. 2(a)). Consequently, 41,616 pixels outside the ROI are reduced to 4,624 pixels. The number of columns of A is reduced to 153 × 153 + 4, 624 = 28, 033. Let the number of detector cells be 255 and the number of views be 110. The multi-resolution system matrix A_M is shown in Fig. 2(b), where the blue dots represent non-zero elements. More complicated system matrices can be constructed in similar ways. To our knowledge, the SVD-based tomographic reconstruction method has not been thoroughly investigated due to the massive size of a solid CT system matrix [13, 23]. Compared to our previous SVD algorithm on 1D PI-segments [12], we refer to our method as SVD-based 2D tomography. In the following sections, we perform a comprehensive evaluation to demonstrate the merits of the algorithm in terms of the system flexibility, fast computation, and accurate reconstruction.

3 Algorithms and computer implementation

3.1 Benchmark algorithms and system configuration

In this paper, the system matrices for the SVD reconstruction are generated using the area integral model (AIM) [30]. We evaluated TSVD, VSVD, and GSVD. The performances of the SVD methods were compared to the filtered backprojection algorithm (FBP [3] with appropriate extrapolation [24, 25] to deal with the data truncation), in which the R-L kernel was applied [26]. The simultaneous algebraic reconstruction technique (SART) [27] and total variation (TV) regularized SART (SART-TV) [28] were also applied. All of the experiments were performed on a high-performance workstation with dual 3.10 GHz Intel Xeon CPUs and 128 GB of memory. SVD of the system matrices was implemented in C with the Intel MKL library, which provides a linear algebra package (LAPACK) that includes the SVD of a matrix. The decomposed results were loaded into Matlab for reconstruction and evaluation. For SART and SART-TV, the projection and backprojection were both implemented by matrix-vector multiplication with A because matrix-vector multiplication is fully optimized inMatlab.

3.2 Data acquisition

3.2.1 Numerical simulation

We employed a modified high-contrast Shepp-Logan phantom and a low-contrast FORBILD head [29] phantom (Fig. 3). The two phantoms are compactly supported in a 50.0 cm×50.0 cm rectangular region. The reconstructed image is a 128×128 mesh grid. The distance from the x-ray source to the iso-center is 53.85 cm, and the distance from the source to the detector is 94.67 cm. An equiangular fan-beam geometry is applied. The fan angle is 0.9592 radians and covers the entire object, and we decrease the fan angle to 0.4796 radians to simulate an interior scan. In our numerical simulations, the standard scan (the entire object was inside the FOV) and interior scan with a circular trajectory are applied using overdetermined system matrices. The number of detector cells and the number of views are both 128.

For the Shepp-Logan phantom, projections were generated by analytically computing overlapping lengths between each x-ray path and ellipsoids. Inconsistencies exist between the analytically calculated intersection lengths and the system matrix coefficients. These inconsistences were viewed as discretization errors. The discretization errors could also be simulated by using different models to generate the projection and SVD-based reconstructions. To further evaluate the discretization errors, the distance-driven model [31 –33] is adopted to generate projections of the FORBILD head phantom. Poisson distributions were used to model the number of photons on each detector cell [34]. Five, 50, and 200 (×10⁴) photons per detector cell were used to simulate air scans for different noise levels without a bowtie filter.

3.2.2 Real physical phantom

A physical phantom is also used to evaluate the SVD methods. The original projection data were acquired on a GE Discovery CT750 HD scanner with a circular scanning trajectory with a radius of 538.5 mm. This geometry defines a FOV with a radius of 249.2 mm. For one rotation, 984 projections are evenly sampled. A total of 888 detector cells is used. The sinogram of the central slice from the equiangular fan-beam geometry is obtained after appropriate pre-processing. The memory size would not allow us to occupy all of the projection data to perform the SVD-based reconstructions. Therefore, the projection data are downsampled. In case 1, the number of detector cells is reduced to 222 by averaging every four adjacent detector cells, and the number of views is downsampled to 246 by discarding three of every four views. In case 2, the number of detector cells is reduced to 444 by averaging every two adjacent detector cells, and the number of views is downsampled to 123 by discarding 7 of every 8 views.

3.3 Quantitative evaluation methods

The root mean square error (RMSE) and structural similarity (SSIM [35]) are applied to quantitatively evaluate the reconstructed images inside the ROI. The RMSE in our experiment is calculated by $RMSE = \sqrt{\sum_{j \in ROI, j \notin Known} {(f_{j}^{r} - f_{j}^{*})}^{2} / N_{R}},$ (23) where $f_{j}^{*}$ is the ground truth, $f_{j}^{r}$ is the reconstructed attenuation coefficient, and N_R is the number of unknown pixels inside the ROI. In addition, the maximum error and standard deviation (STD) inside the rectangular areas marked E and G in Fig. 3 are also measured in several experiments.

4 Results

4.1 Shepp-Logan phantom reconstruction results

4.1.1 Standard scan geometry

The Shepp-Logan phantom is reconstructed from a standard scan dataset. A total of 128 views are uniformly distributed in the scan range [0, 2π]. The parameters are r = 71 % (ɛ = 21.3) for TSVD and ξ = 3.311 for VSVD. The maximum numbers of iterations for SART and SART-TV are 300. The gradient descent method is applied to minimize the TV, and the step sizes are carefully tuned. Reconstructed images are shown in Fig. 4. A narrow display window [0.15, 0.25] is used to determine the differences between the reconstructed images from the proposed and benchmark methods. We observe that FBP reconstructs the worst images with severe artifacts around the object support. SART-TV reconstructs smooth images. SART-TV can visually eliminate artifacts around the object; however, the parameters must be carefully tuned to avoid removing the high frequency information (see the three tiny ellipsoids at the bottom of the Shepp-Logan phantom in Fig. 4). Optimizing the parameters for regularized iterative reconstruction algorithms is a challenging problem. Note that the TV minimization assumes a piecewise constant image model, which is more suitable for reconstruction of the Shepp-Logan phantom. It is not surprising that SART-TV reconstructs images with the smallest maximum error, average error, and STD in E. Both TSVD and VSVD can suppress streak artifacts somewhat compared to FBP, and they retain more high frequency information than SART/SART-TV. TSVD can suppress artifacts outside the object support better than SART, while more high frequency artifacts occur inside the object support than with SART. The VSVD reconstruction has a similar image quality to SART but with more high frequency information. GSVD outperforms TSVD and VSVD. GSVD not only suppresses artifacts around the object support but also retains more high-frequency information with smooth images than TSVD and VSVD. However, we can observe undershoots inside the skull. After evaluating hundreds of values of ξ, undershoots are inevitable under the smallest RMSE/largest SSIM criteria. The RMSE and SSIM are unsatisfying if we vary other parameters to avoid undershoots. The quantitative evaluations (Table 1) show that GSVD provides the smallest RMSE. Although GSVD is inferior to SART-TV in terms of the criteria described above, SART-TV requires 300 iterations to outperform GSVD, which is time-consuming. Although many iterative reconstruction algorithms can converge to the optimal solution of the objective function in tens of iterations, they are still slower than the SVD-based methods, which only require one matrix-vector multiplication, for small image reconstructions.

TSVD is applied to reconstruct images from noise-free projections by keeping different magnitude ratios of the singular values. Figure 5 shows the reconstructed results from keeping r = 10% , 20 % , 30 % , 60 % , 80 % , and 90 % of the singular value magnitudes. A very blurred contour of the Shepp-Logan phantom is produced by keeping 10% of the singular value magnitudes. The overall shape of the Shepp-Logan phantom can be observed with 30% of the singular values, but it contains visible ring artifacts. With more singular values, additional details can be reconstructed; however, artifacts caused by the inconsistencies also become significant. When all of the singular values are involved, artifacts cover all of the image information. We also evaluate VSVD with hundreds of different values of ξ. Figure 6 shows 6 representative results from different values ofξ. With larger values of ξ, the image becomes smoother with less high frequency information, and there is a DC down-shift trend.

4.1.2 Interior scan geometry

The experiments described above are repeated with the interior scan geometry. The fan angle is decreased to 0.4796 radians, and the detector cells and views remain the same. The same-sized Shepp-Logan phantom is reconstructed. Because of data truncation, the projections are preprocessed according to [24] and [25] for FBP with appropriate extrapolation. Similarly, four noise levels (no noise, 5 × 10⁴, 50 × 10⁴, and 200 × 10⁴ photons per detector cell) are applied. The parameters are r = 0.71 (ɛ = 21.3) for TSVD, ξ = 3.4 for VSVD and ξ = 1.9 for GSVD. Figure 7 shows the reconstructed ROIs from the noise-free/noisy datasets in a narrow display window [0.15, 0.25]. The FBP results show little because the extrapolation downshifts the DC. Without regularization, SART provides visually satisfying results with truncation artifacts. In contrast, SART-TV reconstructs smoother images but eliminates higher frequency information when the parameters are not tuned carefully. Two small ellipsoids at the center of the FOV are blurred because of the regularization in SART-TV. However, the noise reduction advantage of SART-TV makes the image piecewise constant and more similar to the ground truth. Both TSVD and VSVD can reconstruct the FOV very well but result in bright truncation artifacts around the FOV edge. This is a result of the inconsistency between the analytical projection calculation and the discrete system matrix estimation. Truncating small singular values in TSVD or modifying reciprocals of singular values in VSVD fluctuates the reconstructed values around the FOV boundary. GSVD outperforms all of the other methods. It faithfully reconstructs the high frequency information (see small ellipsoids at the center of the Shepp-Logan phantom) and suppresses truncation artifacts around the FOV. The quantitative evaluations are summarized in Table 2. With the exception of SART-TV, which has the best STD, GSVD provides the best quantitative image quality of all of the benchmark methods. The constant shift of FBP makes it worst in terms of the RMSE. TSVD slightly outperforms VSVD in terms of RMSE. Both TSVD and VSVD suppress high frequency information but in different ways; TSVD suppresses the high frequency information by abandoning the corresponding singular values, whereas VSVD weakens the high frequency information by modifying the magnitudes of the reciprocals of the singular values.

4.2 FORBILD head phantom reconstruction

The 128×128 FORBILD head phantom is reconstructed to evaluate the SVD methods for low-contrast objects from an interior scan. In many real applications, the ROI is not accurately located at the center of the object. In this case, we relocate the ROI, extrapolate the projection dataset, and enlarge the reconstructed image size to cover the entire object. A standard scan is performed with a 128-element detector, and 271 views are acquired. Poisson noise is added to simulate four different noise levels. We also apply the extrapolation to deal with the projection truncation for FBP [24, 25]. The parameters are r = 90% (ɛ =27.2) for TSVD, ξ = 4.4 for VSVD and ξ=1.5 for GSVD. Figure 8 shows the reconstructed ROIs in a [0.8, 1.5] display window. After extrapolation, FBP can reconstruct a visually acceptable image. All of the methods can reconstruct the ROI except for low contrast eyeball details because of noise and data truncation. GSVD reconstructed the best image quality without truncation artifacts, while the other non-SVD-based methods require extrapolation or regularization to eliminate severe truncation artifacts. Furthermore, only GSVD accurately reconstructs the FORBILD head phantom boundary. The quantitative evaluations of the reconstructed results are summarized in Table 3 and show that GSVD provides the best image quality compared to the other benchmark methods. SART-TV reconstructs the ROI with the smallest STD in region G. When the ROI only occupies a very small region of the object support, the truncation artifacts in SART and SART-TV are significant. The SVD methods can reduce truncation artifacts very well even without prior information or extrapolation. Both TSVD and VSVD reconstruct images with more high frequency information and noise. When the Poisson noise is high, VSVD slightly outperforms TSVD.

4.3 Physical phantom reconstruction

A real physical phantom is reconstructed in a 150×150 grid that covers an area of 500.0 mm×500.0 mm with a pixel size of 3.33 mm×3.33 mm. Because there is no ground truth and this phantom is actually piecewise constant, we apply SART-TV with carefully tuned parameters to reconstruct the standard image from all of the available projections. The reconstructed results from the different methods are shown in Fig. 9. The display window is [–1000, –500] HU. In case 1, although FBP reconstructs the image with visible artifacts, it has the best ability to resolve the details (small tubes inside red circles). Visually, the tiny tube structures can be maintained in TSVD and VSVD at the cost of greater streak artifacts. Both SART and GSVD slightly blur the small tubes. SART has the best performance at suppressing streak artifacts. The SVD methods can also inhibit streak artifacts, but their performance is not as good as SART. Only GSVD is comparable to SART. In case 2, with fewer views, FBP reconstructs the image with more significant artifacts. However, all of the other methods can clearly reconstruct the tiny tube structures, although they are not as obvious as in FBP. We can observe streak artifacts in the piecewise part inside the phantom in the results from SART, TSVD, and VSVD. These streak artifacts are suppressed by the regularization in GSVD. However, the tube structures are slightly blurred in GSVD. The quantitative results in Table 4 show that the smallest RMSE was obtained by SART in case 1. However, the VSVD-reconstructed image has the best SSIM. In case 2, GSVD has the best performance in both RMSE and SSIM.

4.4 SVD methods with a Known Sub-region in the ROI

Equations (17) to (21) show that the interior problem can be solved uniquely and stably by the SVD methods when a sub-region is known inside the ROI. SVD only needs to be performed on A_Ω∖K. Because the compact support of the object is usually known or measurable, we can treat pixels outside the compact support as the known part. In this experiment, a 215×215 image matrix is reconstructed from noise-free projections. The known sub-region is a 64.6 mm×62.5 mm rectangular region (Fig. 3(a)) that includes 24,086 known pixels and 22,139 unknown pixels. There are 324 detector cells, and 60 views are uniformly sampled in a standard scan range.

When there is no discretization error and noise, SVD can exactly reconstruct the image (Fig. 10(a)). The reconstruction error in a much narrower display window is shown in Fig. 10(b). Although the top and bottom parts of the image outside the ROI are noisy, the pixels inside the compact support are accurate. The horizontal striped errors can be viewed as numerical limitations. Figure 10(a) can be regarded as an exact reconstruction. Compared to the SVD methods, SART cannot accurately reconstruct the entire image (Fig. 10(c)). Both the TSVD and VSVD methods are applied for noisy projections. We also solved Equation (7) with the same parameter ξ for VSVD in an iterative fashion, which is referred to as GD-Tik. The aforementioned projections with six noise levels are applied for the evaluation. Threshold values of ɛ (0.063, 0.060, 0.058, 0.0571, 0.0565 and 0.0560 for six noise levels) and the parameter ξ (0.0527, 0.0502, 0.0483, 0.0467, 0.046 and 0.0454 for six noise levels) for VSVD are empirically chosen. A maximum of 300 iterations is used for GD-Tik and SART. The known sub-region is enforced in every iteration. The interior reconstruction results are shown in Fig. 11. A wider display window [0, 0.5] is applied to better show the reconstructed images with the known part inside the FOV. Artifacts around the boundaries of the known sub-region are visible in all of the reconstructed results, and they are more obvious in the results from the iterative methods. They are caused by the high frequency disturbances of the jumps between the known and unknown parts and the sensitivity of human visual perception to image edges. Figure 12 shows the quantitative evaluations. VSVD outperforms TSVD with all noise levels. Because the VSVD results can be viewed as the results of GD-Tik with an infinite number of iterations, the performance of GD-Tik is comparable to that of VSVD in some criteria. Although SART outperforms TSVD and VSVD, for the case of 5 × 10⁴ photons, it requires at least 24 (6.03 seconds) iterations to outperform TSVD and 55 (12.54 seconds) iterations to outperform VSVD. For the case of 200 × 10⁴ photons, SART requires 43 iterations (10.03 seconds) to outperform TSVD and 101 iterations (22.53 seconds) to outperform VSVD. GD-Tik requires 144 (15.82 seconds) and 321 (35.32 seconds) iterations to minimize Equation (6) for 5 × 10⁴ and 200 × 10⁴ photons, respectively. Combining the merits of the SVD methods and SART, we can use results with more high frequency information and noise from the SVD methods to initialize the iterative reconstruction algorithms for the best performance in terms of the image quality and computational cost.

4.5 SVD methods for multiple-resolution images

We evaluate SVD methods with a multi-resolution strategy as described in Equation (22). Because GSVD requires twice the memory to perform and we were limited by the memory of our current workstation, we did not apply GSVD in this experiment. The inconsistencies between the matrices for the projection generation and image reconstruction are treated as discretization errors. In this study, VSVD is used (ξ = 2.98) to calculate the pseudo inverse of A_M. The reconstructed results are shown in Fig. 13. As we expect, the images within the ROI can be accurately reconstructed with high resolution. The pixels outside the ROI are not important for the interior tomography.

Similarly, projections with six noise levels are evaluated in the multiple-resolution TSVD and VSVD. We chose ɛ = 21.1 and ξ = 2.3 for all noise levels. Two multiple-resolution-based reconstruction methods are also evaluated using a known sub-region inside the ROI. The reconstructed results are shown in Fig. 14, and the quantitative evaluation results are listed in Table 5. Although the differences in image quality are not obvious in Fig. 14, the results in Table 5 clearly show that VSVD outperforms TSVD.

5 Discussion and conclusion

In this paper, we proposed and evaluated SVD-based CT reconstruction methods. GSVD outperforms the other benchmark methods in both standard scan and interior scan cases. The SVD methods can accurately reconstruct the ROI regardless of where it is located. An overdetermined system matrix is necessary for an accurate reconstruction. For noise-free projections, if the system matrix is overdetermined and the pixel is illuminated by at least one x-ray path, the pixel can be accurately reconstructed even if it is outside the ROI. For noisy projections, the ROI can be accurately reconstructed by GSVD, TSVD, and VSVD. However, the pixels outside the ROI can no longer be accurately reconstructed. These conclusions are consistent with theoretical results from [6] and [20]. When the matrix was overdetermined, Hansen derived perturbation bounds by establishing the stability of the SVD method to recover all of the variables. When the number of measurements is smaller than the number of unknown variables (pixels), the system matrix is underdetermined, Hansen’s results cannot be directly applied [21], and SVD does not guarantee accurate reconstruction even for noise-free projections. However, if the object support and a sub-region are known, the system matrix can be grouped into known and unknown parts. If the system matrix that corresponds to the unknown part is overdetermined, the SVD methods allow accurate reconstructions of noise-free projections. Meanwhile, the pseudo-inverse of the singular value matrix can be modified to obtain a regularized solution.

Currently, the SVD of a fairly large matrix is still a difficult problem due to the high computational complexity and the large memory requirements. To resolve the interior issue of a high-resolution image, we propose to represent the image using a multi-resolution scheme. This method reduces the size of the system matrix and provides a practical way to reconstruct a high-resolution ROI with limited memory capacity. The proposed methods provide a fast image reconstruction tool. SVD only needs to be performed once when the CT geometry is fixed.

In our experiments, the SVD methods had a lower computational cost than the other methods. Reconstructing a 128×128 image required approximately 0.2 seconds. For SART, it consumed 63.7 seconds for 300 iterations. Furthermore, the computational costs of the projection and backprojection operations depend on the scanning geometry configuration, the projection and backprojection models, and their implementation efficiency. The computational time of the SVD methods was also less than FBP in our experiments because the SVD methods can fully utilize the computational resources to accelerate matrix-vector multiplication. Although the SVD methods have greater computational complexities than FBP, FBP does not have the flexibility to deal with interior scans and other irregular scanning geometries, and the system matrix can be freely modified to perform different reconstruction tasks, including ROI reconstruction.

In conclusion, SVD methods were developed for fast and accurate reconstruction of a broad class of inverse problems that have no analytical solutions. The reconstruction can be rapidly simplified to matrix-vector multiplication. The system matrix can be compressed in both the column and row dimensions by representing the image in a multi-resolution scheme. Additional theoretical analysis and extensive experiments will be conducted to further improve both the SVD performance and the quality of the reconstructed image.

Footnotes

Acknowledgments

This work was partially supported by NSF CAREER Award 1540898 and NSF Collaborative Award 1619550.

References

Linton

O.W.

and Mettler

F.A.

, National Conference on Dose Reduction in CT, with an Emphasis on Pediatric Patients, Am J Roentgenol 181(2) (2003), 321–329.

Wang

and Yu

, The meaning of interior tomography, Phys Med Biol 58(16) (2013), R161.

Natterer

, The Mathematics of Computerized Tomography. Philadelphia: SIAM: Society for Industrial and Applied Mathematics, 2001.

Kudo

, Courdurier

, Noo

and Defrise

, Tiny a priori knowledge solves the interior problem in computed tomography, Phys Med Biol 53(9) (2008), 2207–2231.

Courdurier

, Noo

, Defrise

and Kudo

, Solving The Interior Problem of Computed Tomography Using A Priori Knowledge, Inverse Probl 24(6) (2008), 65001.

, Yu

, Wei

and Wang

, A general local reconstruction approach based on a truncated Hilbert transform, J Biomed Imaging 2007(1) (2007), 2.

, Yu

and Wang

, Exact Interior Reconstruction with Cone-Beam CT, Int J Biomed Imaging 2007 (2007), pp. 1–5.

and Wang

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

Gel’fand

I.M.

and Graev

M.I.

, Crofton’s function and inversion formulas in real integral geometry, Funct Anal Its Appl 25(1) (1991), 1–5.

10.

, Yu

and Wang

, Gel’fand–Graev’s reconstruction formula in the 3D real space, Med Phys 38(Suppl 1) (2011), S69–S75.

11.

, Ye

and Wang

, Interior Reconstruction Using the Truncated Hilbert Transform via Singular Value Decomposition, J X-Ray Sci Technol 16(4) (2008), 243–251.

12.

Jin

, Katsevich

, Yu

, Wang

, Li

and Chen

, Interior Tomography With Continuous Singular Value Decomposition, IEEE Trans Med Imaging 31(11) (2012), 2108–2119.

13.

Sun

, Ma

and Jin

, SVD reconstruction algorithm in 3D SPECT imaging, IEEE Nucl Sci Symp Conf Rec 6 (2004), 3527–3530.

14.

Smith

M.F.

, Floyd

C.R.

, Jaszczak

R.J.

and Coleman

R.E.

, Reconstruction of SPECT images using generalized matrix inverses, IEEE Trans Med Imaging 11(2) (1992), 165–175.

15.

Zeng

G.L.

and Gullberg

G.T.

, An SVD study of truncated transmission data in SPECT, IEEE Trans Nucl Sci 44(1) (1997), 107–111.

16.

Shim

and Cho

, SVD pseudoinversion image reconstruction, IEEE Trans Acoust Speech Signal Process 29(4) (1981), 904–909.

17.

Verhoeven

, Limited-data computed tomography algorithms for the physical sciences, Appl Opt 32(20) (1993), 3736–3754.

18.

VanderPlas

J.T.

, Connolly

A.J.

, Jain

and Jarvis

, Three-dimensional Reconstruction of the Density Field: An SVD Approach to Weak-lensing Tomography, Astrophys J 727(2) (2011), 118.

19.

Bhuyan

, Singh

S.B.

and Bhuyan

P.K.

, Application of generalized singular value decomposition to ionospheric tomography, Ann Geophys 22(10) (2004), 3437–3444.

20.

Hansen

, Truncated Singular Value Decomposition Solutions to Discrete Ill-Posed Problems with Ill-Determined Numerical Rank, SIAM J Sci Stat Comput 11(3) (1990), 503–518.

21.

Hansen

P.C.

, Discrete inverse problems: Insight and algorithms, Philadelphia: Society for Industrial and Applied Mathematics, 2010.

22.

Golub

G.H.

and Loan

C.F.V.

, Matrix Computations, Fourth edition edition. Baltimore: Johns Hopkins University Press, 2012.

23.

Buzug

T.M.

, Computed Tomography: From Photon Statistics to Modern Cone-Beam CT. Springer Science & Business Media, 2008.

24.

Hsieh

, Chao

, Thibault

, Grekowicz

, Horst

, McOlash

and Myers

T.J.

, A novel reconstruction algorithm to extend the CT scan field-of-view, Med Phys 31(9) (2004), 2385–2391.

25.

Sourbelle

, Kachelriess

and Kalender

W.A.

, Reconstruction from truncated projections in CT using adaptive detruncation, Eur Radiol 15(5) (2005), 1008–1014.

26.

Ramachandran

G.N.

and Lakshminarayanan

A.V.

, Three-dimensional Reconstruction from Radiographs and Electron Micrographs: Application of Convolutions instead of Fourier Transforms, Proc Natl Acad Sci 68(9) (1971), 2236–2240.

27.

Wang

and Jiang

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

28.

, Ji

and Wang

, SART-Type Image Reconstruction from Overlapped Projections, Int J Biomed Imaging 2011 (2011), 1–7.

29.

“http://www.imp.uni-erlangen.de/phantoms.”

30.

and Wang

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

31.

De Man

and Basu

, Distance-driven projection and backprojection in three dimensions, Phys Med Biol 49(11) (2004), 2463–2475.

32.

Basu

and De Man

, Branchless distance driven projection and backprojection, 2006, pp. 60650Y–60650Y–8.

33.

Liu

, Fu

, De Man

and Yu

, GPU Acceleration of Branchless Distance Driven Projection and Backprojection, presented at the 4th International Conference on Image Formation in X-Ray Computed Tomography, Bamberg, Germany, 2016.

34.

Hsieh

, Computed tomography: Principles, design, artifacts, and recent advances 2rd edition, 114. Society of Photo-Optical, 2009.

35.

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.