Ultra-limited-angle CT image reconstruction algorithm based on reweighting and edge-preserving

Abstract

BACKGROUND:

Ultra-limited-angle image reconstruction problem with a limited-angle scanning range less than or equal to $\frac{π}{2}$ is severely ill-posed. Due to the considerably large condition number of a linear system for image reconstruction, it is extremely challenging to generate a valid reconstructed image by traditional iterative reconstruction algorithms.

OBJECTIVE:

To develop and test a valid ultra-limited-angle CT image reconstruction algorithm.

METHODS:

We propose a new optimized reconstruction model and Reweighted Alternating Edge-preserving Diffusion and Smoothing algorithm in which a reweighted method of improving the condition number is incorporated into the idea of AEDS image reconstruction algorithm. The AEDS algorithm utilizes the property of image sparsity to improve partially the results. In experiments, the different algorithms (the Pre-Landweber, AEDS algorithms and our algorithm) are used to reconstruct the Shepp-Logan phantom from the simulated projection data with noises and the flat object with a large ratio between length and width from the real projection data. PSNR and SSIM are used as the quantitative indices to evaluate quality of reconstructed images.

RESULTS:

Experiment results showed that for simulated projection data, our algorithm improves PSNR and SSIM from 22.46db to 39.38db and from 0.71 to 0.96, respectively. For real projection data, our algorithm yields the highest PSNR and SSIM of 30.89db and 0.88, which obtains a valid reconstructed result.

CONCLUSIONS:

Our algorithm successfully combines the merits of several image processing and reconstruction algorithms. Thus, our new algorithm outperforms significantly other two algorithms and is valid for ultra-limited-angle CT image reconstruction.

Keywords

Computed Tomography (CT)ultra-limited-angle reconstruction condition number

1 Introduction

In Computed Tomography (CT), the limited-angle image reconstruction is to invert the two-dimensional image function f (x, y) from the following limited-angle Radon transform [1], $Rf (p, θ) = \int_{ℓ : p = x cos θ + y sin θ} f (x, y) d ℓ, 0 \leq θ \leq θ_{0} < π$ (1) where p is the distance the X-ray ℓ from the origin on the plane R², θ is the counterclockwise angle of the normal vector of the ray and the horizontal axis of plane R², and Rf (p, θ) is the given projection data.

Since the 1980s, some progress of the limited-angle CT image reconstruction problem has been achieved. However, it is still very tough to achieve a valid reconstruction from a scanning angle range less than or equal to $\frac{π}{2}$ . In [1, 2], Smith et al. studied the uniqueness and stability of the solution to the limited-angle CT image reconstruction problem. In [3], Louis explained theoretically that the ill-posedness of limited-angle problem is aggravated with the increase of missed angle and showed numerically that some parts of image spectrum cannot be irrecoverable when the projection angle range is less than or equal to $\frac{π}{2}$ . In [4], Frikel et al. explained why the streak artifacts are present at the end of the limited-angle range by using micro-local analysis. In [5], Quinto indicated that the image boundaries tangent to the projection lines in the limited data set should be easy to reconstruct, such boundaries are called as visible boundaries, otherwise, the image boundaries should be difficult to reconstruct and are called as invisible boundaries. Traditional image reconstruction algorithms require complete projection data. For the limited-angle projection data, the existing image reconstruction algorithms include mainly the Total Variation (TV)-based algorithms and deep learning. In [6], Sidky et al. applied the TV algorithm to the incomplete data CT image reconstruction problems. In [7], Yu et al. designed simultaneous algebraic reconstruction technique with TV regularization (SART-TV), In [8], Chen et al. presented SART-ATV reconstruction approach based on anisotropic TV (ATV) minimization method. In [9], Xu et al. gave the alternating edge-preserving diffusion and smoothing (AEDS) algorithm by introducing the regularization term along the horizontal-direction and the regularization term vertical-direction, respectively. These TV-based methods first got an intermediate result by traditional iterative reconstruction algorithms (such as ART and SART), and then recovered invisible boundaries based on visible boundary information of the intermediate result.

However, the fundamental difficulty caused by the considerably large condition number remains. In [10 –13], several different deep learning architectures were designed for the limited-angle CT image reconstruction problem. However, these deep learning technologies always need tremendous training data to predict the weight factors of the reconstruction architectures. These algorithms are mainly used for the limited-angle reconstruction problem with a limited-angle scanning range more than $\frac{π}{2}$ . In the industrial and medical fields, due to the geometric limitation of scanning system and the dose requirement of radiation, projection data can only be collected in a fairly small angle range. For example, digital breast tomography reduces the ray dose mainly by obtaining less projection data from a limited-angle scanning range less than or equal to $\frac{π}{2}$ [14]. The limited-angle CT image reconstruction with a scanning range less than or equal to $\frac{π}{2}$ is called as the ultra-limited-angle CT image reconstruction problem [15]. The ultra-limited-angle CT image reconstruction problem is severely ill-posed. In [15], Li et al. designed the Sinogram-Inpainting-GAN (SI-GAN) network architecture to restore missing projection data and developed SI-GAN + SART-TV algorithm for the ultra-limited-angle CT image reconstruction. However, as with other deep learning technologies, SI-GAN + SART-TV algorithm also needs tremendous training data. In [16], Zhang et al. formulated the limited-angle problem as a convex optimization program with directional TV, and developed an iterative algorithm to solve the program. The developed algorithm was demonstrated to reconstruct accurately images from data over a significantly reduced angular range. At present, the ultra-limited-angle CT image reconstruction problem is of importance and challenge.

The objective of this paper is to get a valid algorithm for the ultra-limited-angle reconstruction problem. The reconstruction equation is the linear system derived by discretizing the limited-angle Radon transform (1) with $θ \in [0, \frac{π}{2}]$ . The theoretical severe ill-posedness of problem leads the linear system to be severely ill-posed and has a considerably large condition number [3]. With the increase of missed angle, the ill-posedness of problem is aggravated and the condition number becomes larger. The difficulties of ultra-limited-angle image reconstruction problem mainly come from two aspects. One is that the considerably large condition number makes the iteration images contain heavy artifacts. The other one is that the loss of large-range angle projection data leads to a lot of invisible boundaries in the reconstructed images. In [17], we gave a reweighted method to improve the condition number of the coefficient matrix of linear system and the corresponding Preconditioned Landweber iteration scheme (Pre-Landweber), and the Pre-Landweber scheme has also been applied elementarily for the limited-angle reconstruction from simulated projection data. In [9], based on the property of image sparsity, the AEDS algorithm recovered partially invisible boundaries and reduce striping artifact from the image by edge-preserving diffusion and smoothing processing. However, the AEDS algorithm does not improve essentially the large condition number of reconstruction equation. Hence, the ability to improve the results is very limited. In this paper, we incorporate the reweighted method into the idea of AEDS algorithm, and then propose an optimized reconstruction model and Reweighted AEDS (Rw-AEDS) algorithm. Each iteration of Rw-AEDS algorithm consists of three steps: step 1 is to get a subsolution by using the Pre-Landweber algorithm and non-negative processing, and steps 2 and 3 perform regularization processing in the horizontal and vertical direction of the subsolution, respectively. In the step 1, the large condition number are improved greatly, and then the subsolution has less artifacts and more visible boundaries which is very important for the next regularization processing.

In contrast with those existing limited-angle image reconstruction algorithms, our algorithm not only improves greatly the considerably large condition number but also uses the prior information to further improve the reconstructed result. Thus, reconstructed results are improved. Based on the simulated projection data with signal to-noise ratio (SNR) of not more than 30db and real projection data, the reconstruction experiments were implemented. The quantitative indices Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity (SSIM) are used to evaluate the quality of reconstructed images. The experimental results show that for the simulated projection data, our algorithm improves PSNR and SSIM from 22.4603db to 39.3844db and from 0.7081 to 0.9606, respectively; for real projection data, our algorithm yields the highest PSNR and SSIM of 30.8916db and 0.8764, and gets a valid reconstructed result. Hence, the proposed Rw-AEDS algorithm is efficient and outperforms significantly the AEDS and Pre-Landweber algorithms in the ultra-limited-angle CT image reconstruction problem.

The rest of this paper is organized as follows. In Section 2, we give some preliminary results. In Section 3, we give a new algorithm for ultra-limited-angle CT image reconstruction. In Section 4, we implement reconstruction experiments. In Section 5, we give a discussion. In Section 6, we conclude this paper.

2 Preliminaries

The continuous various image f (x, y) has compact support in [0, 1] × [0, 1]. By Partitioning [0, 1] × [0, 1] into N = N₁ × N₁ square cells on an average, the continuous image function f (x, y) can be discretized approximately as f (x_i, y_j), where i, j = 1, ⋯ , N₁. Let k = (i - 1) N₁ + j, f (x_i, y_j) can be rearranged as a one-dimensional vector $x = {(x_{1}, \dots, x_{N})}^{T}$ . Let a_lk be the secant length of the l th ray within the k th cell, where l = 1, ⋯ , M and k = 1, ⋯ , N. Then the limited-angle Radon transform (1) may be approximately as a following linear system $Ax = b$ (2) where A = [a_lk] _M×N and b = (b₁, ⋯ , b_M) ^T. In theory, the limited-angle Radon transform has a unique solution [3]. With the increase of grid number, the discrete one-dimensional vector x and the coefficient matrix A tend to the continuous image function f (x, y) and the limited-angle Radon transform (1) in a particular form, respectively. If the grid is enough, we assume that the linear system (2) has the unique solution $\bar{x}$ . Let A ₀ = A ^T A and b ₀ = A ^T b , the normal equation of linear system (2) is $A_{0} x = b_{0},$ (3) where the matrix A ₀ is symmetric positive definite. Let ${μ_{0, k}}_{k = 1}^{N}$ with the algebraic multiplicity p_k and ${u_{j, k}}_{k = 1}^{N} (1 \leq j \leq p_{k})$ be the distinct eigenvalues of A ₀ and the corresponding orthonormal eigenvectors, respectively. These eigenvalues satisfy μ_0,1 > μ_0,2 > ⋯ > μ_0,N > 0. Let I _pk be the pk-order unit matrix and U = ( U ₁, ⋯ , U _N), where U _k = ( u _{1
,k}, ⋯ , u _{pk
,
k}). Then matrix A ₀ can be diagonalized into $U^{T} A_{0} U = Diag (μ_{0, 1} I_{p 1}, μ_{0, 2} I s_{p 2}, \dots, μ_{0, N} I_{pN}) .$ (4)

Definition 2.1 For a positive definite matrix B, let μ_max and μ_min be its maximum eigenvalue and minimum eigenvalue, respectively. Then, condition number of B is defined as $κ (B) : = \frac{μ_{max}}{μ_{min}}$ .

By introducing the two regularization terms along the horizontal-direction and vertical-direction respectively, a reconstruction model is given and as follows [9] $x = arg min {∥ Ax - b ∥_{W}^{2} + λ_{1} φ (E_{n} x) + λ_{2} ψ (S_{v} x)}$ (5) where $∥ Ax - b ∥_{W}^{2}$ is the fidelity term, φ (E_n x ) and ψ (S_v x ) are regularization terms along the horizontal-direction and vertical-direction, respectively. The forms of φ (E_n x ) and ψ (S_v x ) are as follows, ${\begin{matrix} φ (E_{n} x) = {∥ D_{n} x ∥}_{0} ≜ \sum_{i = 1}^{N_{1}} sgn (| f_{j, i} - f_{j, i - 1} |) \\ ψ (S_{v} x) = {∥ S_{v} x ∥}_{0} ≜ \sum_{j = 1}^{N_{1}} sgn (| f_{j, i} - f_{j - 1, i} |) \end{matrix}$ (6) where sgn () is the sign function such that sgn (0) =0. AEDS algorithm is an alternative iterative approach for solving the reconstruction model (5) [9].

AEDS algorithm
Input: Initial guess x ⁽⁰⁾, projection data b , maximum iteration number
s1, parameter λ₁, λ₂, relaxation coefficient ω_s = 0.5
Output: x ^(s)
for s = 1:s1 do
Step 1:
$x^{(s + \frac{1}{3})} = SART (x^{(s)}, b, ω_{s})$ , $x^{(s + \frac{1}{3})} = shrink (x^{(s + \frac{1}{3})}, 0)$
Step 2: Solving optimization problem:
$x^{(s + \frac{2}{3})} = \arg \min {{‖ x^{(s + \frac{1}{3})} - x ‖}_{2}^{z} + λ_{1} ∥ D_{n} x ∥_{0}}$
Step 3: Solving optimization problem:
$x^{(s + 1)} = \arg \min {{‖ x^{(s + \frac{2}{3})} - x ‖}_{2}^{z} + λ_{2} ∥ S_{v} x ∥_{0}}$
end

In above AEDS algorithm, each iteration contains three steps. Step 1 is to get the subsolution x ^(s+1/3) by using SART algorithm and non-negative processing. Steps 2 and 3 perform regularization processing in the horizontal and vertical direction of the subsolution x ^(s+1/3), respectively. Thus, the step 1 is very critical for AEDS algorithm. However, for the ultra-limited-angle image reconstruction problem with a limited-angle scanning range less than or equal to $\frac{π}{2}$ , the considerably large condition number κ ( A ₀) makes the subsolution x ^(s+1/3) obtained by SART algorithm have numerous invisible boundaries, and then affects the final reconstruction result of AEDS algorithm. Hence, it is necessary to improve the considerably large condition number κ ( A ₀).

3 Methods

3.1 Reweighted method

By the orthonormal transformation \bm U in equality (4), for $a_{0} \geq 1 + \frac{1}{κ (A_{0})}$ , then $U^{T} (a_{0} μ_{0, 1} I - A_{0}) U = Diag ((a_{0} μ_{0, 1} - μ_{0, 1}) I_{p 1}, \dots, (a_{0} μ_{0, 1} - μ_{0, N}) I_{pN}) .$ (7)

From μ_0,1 > μ_0,2 > ⋯ > μ_0,N > 0 and it follows that 0 < a₀μ_0,1 - μ_0,1 < a₀μ_0,1 - μ_0,2 < a₀μ_0,1 - μ_0,N. In order to improve the considerably large condition number κ ( A ₀), by left-multiplying the normal equation (3) by a positive definite matrix $V_{1}^{- 1} = a_{0} μ_{0, 1} I - A_{0}$ and repeating this process for many times, we give the following reweighted method [17], ${\begin{matrix} A_{n} x = b_{n} \\ A_{n} = V_{n}^{- 1} A_{n - 1}, b_{n} = V_{n}^{- 1} b_{n - 1}, V_{n}^{- 1} = a_{n - 1} μ_{n - 1, 1} I - A_{n - 1} \end{matrix}$ (8) where a_n-1 are weighted parameters, μ_n-1,1 are the largest eigenvalues of the matrices A_n-1 and the weighting times n = 1, 2, ⋯.

Theorem 3.1 [17]Assume that $a_{n} \geq 1 + \frac{1}{κ (A_{n})}$ and {a_n} is a convergent sequence with $lim_{n \to \infty} a_{n} \geq 2$ Then κ ( A _n+1) < κ ( A _n) and lim _n→∞ κ ( A _n) =1.

With the reweighted method (8), we get the following Pre-Landweber scheme with the preconditioner ${\bar{V}}_{n} = \prod_{k = 1}^{n} V_{k}^{- 1}$ for solving the linear system (2) [17], $x^{(s + 1)} = x^{(s)} + r_{s} {\bar{V}}_{n} A^{T} (b - A x^{(s)}),$ (9) where s is the iteration number and relaxation coefficient $r_{s} = \frac{2}{1.00002 μ_{n, 1}}$ . Note that ${\bar{V}}_{n} A^{T} b = b_{n}$ and ${\bar{V}}_{n} A^{T} A = A_{n}$ , which can be derived easily by the equality (8).

3.2 Optimization model and Rw-AEDS algorithm

With the reweighted method (8), the reconstruction model (5) may be optimized as $x = arg min {{∥ A_{n} x - b_{n} ∥}_{W}^{2} + λ_{1} φ (E_{n} x) + λ_{2} ψ (S_{v} x)} .$ (10)

Combining the reweighted method (8) of improving condition number and AEDS algorithm, we give the following Rw-AEDS algorithm:

Rw-AEDS algorithm
Input: Initial guess x ⁽⁰⁾, projection data b , maximum iteration number
s1, parameters λ₁, λ₂, relaxation coefficient $r_{s} = \frac{2}{1.00002 μ_{n, 1}}$
Output: x ^(s)
for s = 1 : s1 do
Step 1:
$x^{(s + \frac{1}{3})} =$ Pre-Landweber $(x^{(s)}, b, r_{s}), x^{(s + \frac{1}{3})} = shrink (x^{(s + \frac{1}{3})}, 0);$
Step 2: Solving optimization problem:
$x^{(s + \frac{2}{3})} = arg min {{∥ x^{(s + \frac{1}{3})} - x ∥}_{2}^{z} + λ_{1} ∥ D_{h} x ∥_{0}};$
Step 3: Solving optimization problem:
$x^{(s + 1)} = arg min {{∥ x^{(s + \frac{2}{3})} - x ∥}_{2}^{z} + λ_{2} ∥ S_{v} x ∥_{0}};$
end

In above Rw-AEDS algorithm, step 1 uses the Pre-Landweber scheme (9). Compared to the AEDS algorithm, the condition number κ ( A _n) is much smaller than κ ( A ₀). Thus, the subsolution x ^(s+1/3) of step 1 in the Rw-AEDS algorithm is better than that of step 1 in the AEDS algorithm. The subsolution x ^(s+1/3) of step 1 has the great effect on the final reconstructed result, and steps 2 and 3 are further processing of the subsolution x ^(s+1/3). Hence, the reconstructed result is also improved.

3.3 Weighting times

In the reweighting process of Rw-AEDS algorithm, the matrix products produce the rounding errors and these rounding errors can be accumulated with the increase of weighting times n. Then in the reweighted method (8), the linear system A _n x = b_n becomes $A_{n} (x + δ x) = b_{n} + δ b_{n},$ (11) where δb _n is the accumulated rounding error. By the general theory of the solution for the linear system [18], then $\frac{∥ δ x ∥}{∥ x ∥} \leq κ (A_{n}) \frac{∥ δ b_{n} ∥}{∥ b_{n} ∥} .$ (12)

Obviously, the smaller the value of $κ (A_{n}) \frac{∥ δ b_{n} ∥}{∥ b_{n} ∥}$ , the the higher is the solution accuracy. With the weighting times n increasing, the condition number κ ( A _n) decreases. However, the accumulated rounding error δ b_n becomes larger, then the value of $\frac{∥ δ b_{n} ∥}{∥ b_{n} ∥}$ becomes larger. Therefore, the weighting times n cannot be very large to ensure the validity of the solution, namely a sufficiently small $κ (A_{n}) \frac{∥ δ b_{n} ∥}{∥ b_{n} ∥}$ . In the reconstruction experiments of Section 4, the maximum weighting times n = 5.

3.4 Selection of the weighted parameter and calculation of the maximum eigenvalue

For the proposed Rw-AEDS algorithm, we now analyze the effect of the weighted parameter a_n on the improvement of condition number. For n = 0, 1⋯, the eigenvalues of the matrix A _n are denoted as ${μ_{n, k}}_{k = 1}^{N}$ satisfying μ_n,1 > μ_n,2 > ⋯ μ_n,N > 0. From the expressions (7) and (8), we have $μ_{n + 1, k} = (a_{n} μ_{n, 1} - t) t = \frac{a_{n}^{2} μ_{n, 1}^{2}}{4} - {(t - \frac{a_{n} μ_{n, 1}}{2})}^{2}, where t \in {μ_{n, k}}_{k = 1}^{N} .$ (13)

For $1 + \frac{1}{κ (A_{n})} \leq a_{n} \leq 2, \frac{μ_{n, 1} + μ_{n, N}}{2} \leq \frac{a_{n} μ_{n, 1}}{2} \leq μ_{n, 1}$ , then $μ_{n + 1, 1} \leq \frac{a_{n}^{2} μ_{n, 1}^{2}}{4}$ and μ_n+1,N = (a_nμ_n,1 - μ_n,N) μ_n,N, so $κ (A_{n + 1}) = \frac{μ_{n + 1, 1}}{μ_{n + 1, N}} \leq \frac{a_{n}^{2} μ_{n, 1}^{2}}{4 (a_{n} μ_{n, 1} - μ_{n, N}) μ_{n, N}} = \frac{a_{n}^{2} κ (A_{n})}{4 [a_{n} κ (A_{n}) - 1]} κ (A_{n}) .$ (14)

Since ${[\frac{a_{n}^{2}}{a_{n} κ (A_{n}) - 1}]}^{'} = [\frac{a_{n}}{{[a_{n} κ (A_{n}) - 1]}^{2}}] > 0, \frac{a_{n}^{2} κ (A_{n})}{4 [a_{n} κ (A_{n}) - 1]} \leq \frac{κ (A_{n})}{2 κ (A_{n}) - 1}$ , and $κ (A_{n + 1}) \leq \frac{κ (A_{n})}{2 κ (A_{n}) - 1} κ (A_{n}) .$ (15)

For $a_{n} > 2, \frac{a_{n} μ_{n, 1}}{2} > μ_{n, 1}$ , then μ_n+1,1 = (a_nμ_n,1 - μ_n,1) μ_n,1 and μ_n+1,N = (a_nμ_n,1 - μ_n,N), so $κ (A_{n + 1}) = \frac{μ_{n + 1, 1}}{μ_{n + 1, N}} = [1 - \frac{κ (A_{n}) - 1}{a_{n} κ (A_{n}) - 1}] κ (A_{n}) > \frac{κ (A_{n})}{2 κ (A_{n}) - 1} κ (A_{n}) .$ (16)

From the inequalities (15) and (16), we have $κ (A_{n + 1}) |_{a_{n} \in [1 + \frac{1}{κ (A_{n})}, 2]} \leq \frac{κ (A_{n})}{2 κ (A_{n}) - 1} κ (A_{n}) < κ (A_{n + 1}) |_{a_{n} \in (2, + \infty]} .$ (17)

Hence, the condition number can be improved faster when the weighted parameter a_n is selected on the interval $[1 + \frac{1}{κ (A_{n})}, 2]$ . For the ultra-limited-angle image reconstruction problem, the condition number is very large. With the finite weighting times, the considerably large condition number is improved but still large. In the experiments of Section 4, the weighted parameters a_k = 1.0005 for k = 0, 1 ⋯ , 4.

Next, we consider the calculation of the maximum eigenvalue μ_n,1. From the equation (13) and it follows that $μ_{n, 1} \leq \frac{a_{n - 1}^{2} μ_{n - 1, 1}^{2}}{4} \leq {(μ_{0, 1})}^{2^{n}} \prod_{i = 0}^{n - 1} {(\frac{a_{i}}{2})}^{2^{n - i}}$ . Since the matrix A is highly sparse, we use the power iteration method to calculate quickly an approximation of μ_0,1 and denoted by ${\bar{μ}}_{0, 1}$ such that $| {\bar{μ}}_{0, 1} - μ_{0, 1} | \leq δ$ , where δ is very small. Then $μ_{0, 1} \leq {\bar{μ}}_{0, 1} + δ$ . Let $a_{n} \geq (1 + \frac{1}{κ (A_{n})}) (\frac{{\bar{μ}}_{0, 1} + δ}{{\bar{μ}}_{0, 1}}) 2^{n}$ and this yields $a_{n} {({\bar{μ}}_{0, 1})}^{2^{n}} \prod_{i = 0}^{n - 1} {(\frac{a_{i}}{2})}^{2^{n - i}} \geq (1 + \frac{1}{κ (A_{n})}) {(μ_{0, 1})}^{2^{n}} \prod_{i = 0}^{n - 1} (\frac{a_{i}}{2}) 2^{n - 1} \geq (1 + \frac{1}{κ (A_{n})}) μ_{n, 1} = μ_{n, N} + μ_{n, 1} .$ (18)

The assumption of Theorem 3.1, $a_{n} \geq 1 + \frac{1}{κ ({sA}_{n})}$ , is equivalent to a_nμ_n,1 ≥ μ_n,N + μ_n,1. Thus, the conclusion of improving condition number still holds if replacing μ_n,1 by ${({\bar{μ}}_{0, 1})}^{2^{n}} \prod_{i = 0}^{n - 1} {(\frac{a_{i}}{2})}^{2^{n - i}}$ . In the reconstruction experiments of Section 4, we use ${({\bar{μ}}_{0, 1})}^{2^{n}} \prod_{i = 0}^{n - 1} {(\frac{a_{i}}{2})}^{2^{n - i}}$ instead of μ_n,1 to avoid the large computation of the maximum eigenvalues.

4 Experimental results

In this section, we investigate the performance of the proposed Rw-AEDS algorithm using simulation and real projection data. The simulation projection data was generated by the Radon transform of the original images. In experiments, we demonstrate the effectiveness and efficiency of our algorithm by comparing it with the Pre-Landweber [17] and AEDS [9] algorithms. The maximum weighting times is 5 and a_k = 1.0005 for l = 0, 1 ⋯ , 4. All the experiments were performed with MATLAB on the PC with a 3.4 GHz Intel Core i7 processor and 64 GBs of memory. The software used is MATLAB. The formulas of PSNR and SSIM are as follows [19] $PSNR (f, f^{'}) = 20 {log}_{10} \frac{255}{\sqrt{Σ_{i = 1}^{N_{1}} Σ_{j = 1}^{N_{1}} {(f_{i, j} - f'_{i, j})}^{2} / (N_{1} \times N_{1})}}$ (19) $SSIM (f, f^{'}) = \frac{(2 uf u_{f^{'}} + 0.01) (2 σ f, f^{'})}{(u f^{2} + uf'^{2} + 0.01) (σ f^{2} + σ f'^{2} + 0.03},$ (20) where f and f′ are the reconstructed image and actual image. f_i,j and f′_i,j represent the pixel values of f and f′, respectively. u_f and u_f′ are the mean value of f and f′, σ_f σ_f′ and are the variance of f and f′, and σ_f,f′ is the covariance of f and f′.

4.1 Experiments on simulation data

In this subsection, to reveal the sensitivity of the proposed Rw-AEDS algorithm to noises, we add the White Gaussian noises of different levels to the simulation projection data. The original image is Shepp-Logan phantom with size of 128×128 (as shown in Fig. 1). Parameters λ₁ = 0.00016 and λ₂ = 0.0002. The maximum iteration number is 200. The reconstruction time for 200 iterations is about 1200s. Figure 1 shows the reconstructed results of Rw-AEDS algorithm from the full-angle projection data. Table 1 gives the PSNR and SSIM of reconstructed images from the full-angle projection data. Figure 2 shows the reconstructed results of different algorithms from the limited-angle projection data, where the known angle range is from $\frac{5 π}{18}$ to $\frac{43 π}{60}$ with $\frac{π}{180}$ increments and 381 rays per projection. Figure 3 plots the PSNR curves. Figure 4 shows the comparison of the 64^th row in different images. Table 2 compares the PSNR and SSIM of reconstructed images from the limited-angle projection data. It can be seen from Fig. 2 that the reconstructed images of the proposed Rw-AEDS algorithm are significantly better than those of other two algorithms, and only the Rw-AEDS algorithm can reconstruct the valid images. It can be found from Table 2 that the PSNR and SSIM of reconstructed images become higher as the increase of the SNR of projection data, and the PSNR and SSIM of Rw-AEDS algorithm are significantly higher than the other two algorithms.

Fig.1

The reconstructed results of Rw-AEDS algorithm from the simulated full-angle projection data. The 1st is the original image. The SNR of projection data in the 2nd, 3rd and 4th images are 5db, 15db and 25db, respectively. The weighting times n = 5.

Table 1

PSNR and SSIM of the reconstructed images of Rw-AEDS algorithm from the simulated full-angle projection data

SNR/db	PSNR/db	SSIM
5	27.2756	0.4104
15	40.4710	0.7877
25	49.8021	0.9638

Fig.2

The reconstructed results from the simulated projection data. The 1st column is the SART algorithm, the 2nd column is the AEDS algorithm, the 3rd column is the Pre-Landweber algorithm and the 4th column is the Rw-AEDS algorithm. The SNR of projection data in the 1st and 2nd rows are 25db and 30db, respectively. The known projection angle range is $\frac{5 π}{18} \sim \frac{43 π}{60}$ . The weighting times n = 5.

Fig.3

PSNR of reconstruction results of Shepp– Logan phantom with White Gaussian noises at each iteration.

Fig.4

Horizontal profiles through the middle of the reconstructed images.

Table 2

PSNR and SSIM of the reconstructed images of different algorithms from the simulated limited-angle projection data

Algorithm	SNR = 25db		SNR = 30db
	PSNR/db	SSIM	PSNR/db	SSIM
SART	19.9336	0.6118	19.9478	0.6357
AEDS	22.6708	0.7875	23.0637	0.8013
Pre-Landweber	21.5839	0.6325	22.4603	0.7081
Rw-AEDS	34.9606	0.8762	39.3844	0.9606

4.2 Experiments on real data

In this subsection, we test the performance of the proposed Rw-AEDS algorithm in practical application through a set of real data. Real image is a flat object with a large ratio between length and width (see Fig. 5). This large ratio increases the difficulty of reconstruction since there are very few visible boundaries [9]. This set of real data was provided by the industrial CT system of the Center for Detection and Imaging Engineering of Capital Normal University, where the x-ray source is YXLON-FXE-225.48 and the flat panel detector is Varian PS2520V. Table 3 lists detailed parameters and other information on this set of real data. In [9], this set of real data has been used to test the performance of AEDS algorithm in the limited-angle CT image reconstruction, where the projection angle range is $[\frac{π}{6}, \frac{5 π}{6}]$ with $\frac{π}{360}$ increments and 1920 rays per projection.

Fig.5

Photograph of the real object used for the experiment.

Table 3

System and geometrical scanning parameters

Parameter	Value
Tube voltage	140 kVp
Current	160 mAs
Scanning angle	[0, π]
Scanning angular interval	π/360
Detector unit width	0.127 mm
Number of detector units	1920
X-ray source to the rotation center distance	311.49 mm
X-ray source to detector distance	697.9082 mm

In the following reconstruction experiments, the full-angle reconstruction result obtained by SART algorithm is used as reference image (see the 1st row in Fig. 6). The limited-angle projection data was extracted from the full angle projection data, where the known projection angle range is $[\frac{π}{4}, \frac{67 π}{90}]$ with $\frac{π}{180}$ increments and 1920 rays per projection. The reconstructed images are size of 256×256. However, for better display, they are clipped to the size of 256×100. Parameters λ₁ = 0.00012 and λ₂ = 0.0002. The maximum iteration number is 100. The reconstruction time for 100 iterations is about 600s. Figure 6 shows the reconstructed results from the different algorithms. Table 4 gives the PSNR and SSIM of reconstructed images. It shows that our algorithm can reconstruct clearer structures and sharper edges than other two algorithms. The reconstructed image of Rw-AEDS algorithm contains more details as the yellow region marked in the 1st image of Fig. 6.

Fig.6

The reconstructed results from the real projection data. The 1st row is the Full-angle SART reconstruction. The 2nd, 3rd and 4th rows are the limited-angle reconstructed images of different algorithms: Pre-Landweber, AEDS and Rw-AEDS, respectively. The known projection angle range is $\frac{π}{4} \sim \frac{67 π}{90}$ . The weighting times n = 5.

Table 4

PSNR and SSIM of the reconstructed images from real projection data

Algorithm	PSNR/db	SSIM
Pre-Landweber	21. 3032	0. 6774
AEDS	19. 8891	0. 4758
Rw-AEDS	30. 8916	0. 8764

5 Discussion

The limited-angle image reconstruction problem with a scanning angle range less than π is ill-posed and this ill-posedness becomes aggravated as the increase of missing angle. As a special case of limited-angle problem, the ultra-limited-angle image reconstruction problem with a scanning range less than or equal to $\frac{π}{2}$ is severely ill-posed. The difficulties mainly come from two aspects. One is that the considerably large condition number makes the iteration images contain heavy artifacts. The other one is that the loss of large-range angle projection data leads to a lot of invisible boundaries in the reconstructed images. Hence, to obtain a valid reconstructed result, reconstruction algorithms need to have the advantages of both improving the condition number and exploiting the prior information to improve the reconstructed result.

The Pre-Landweber and AEDS algorithms are two latest limited-angle CT image reconstruction algorithms. They improve the reconstructed results in different ways, where Pre-Landweber algorithm improves the condition number and AEDS algorithm utilizes the visible boundaries. However, for the ultra-limited-angle image reconstruction problem studied in this paper, the reconstructed results of these two algorithms have a low PSNR value (see Tables 2 and 4) and are invalid (see Figs. 2 and 6). Since the difficulties of ultra-limited-angle image reconstruction problem come from two aspects, improving only one of them is not enough to achieve a valid reconstruction. To this end, we propose the Rw-AEDS algorithm by incorporating the reweighted method of improving condition number into the idea of AEDS algorithm. In contrast with those existing limited-angle image reconstruction algorithms, our algorithm not only improves the considerably large condition number but also uses the prior information to further improve the reconstructed result. By implementing the reconstruction from the simulated projection data with SNR of not more than 30db and real projection data, our algorithm is compared with the Pre-Landweber and AEDS algorithms. As shown in Figs. 2 and 6, the reconstructed images of our algorithm contain less artifacts and more details than those of other algorithms, and only our algorithm is valid for ultra-limited-angle image reconstruction problem. From Table 2, it can be found that our algorithm improves PSNR and SSIM from 22.4603db to 39.3844db and from 0.7081 to 0.9606, respectively. From Table 4, our algorithm yields the highest PSNR and SSIM of 30.8916db and 0.8764 for real data. Thus, our algorithm can improve significantly PSNR and SSIM of reconstructed images and is valid for ultra-limited-angle image reconstruction problem.

6 Conclusion

In conclusion, this paper has proposed a new algorithm referred to as the Rw-AEDS algorithm for the ultra-limited-angle CT image reconstruction. This algorithm combines the merits of AEDS and the reweighted method, resulting in much better performance. The AEDS algorithm utilizes the image sparsity, and the reweighted method improves the large condition number of reconstruction equation. The experimental results verify the superiority and effectiveness of the proposed Rw-AEDS algorithm. It is worth noting that our algorithm can also be extended to fan-beam case.

Footnotes

Acknowledgments

This work is partially supported by the National Natural Science Foundations of China (61931003, 61671004). We thank sincerely the Center for Detection and Imaging Engineering of Capital Normal University for providing the experimental data.

References

Smith

and Solmon

, et al., Practical and mathematical aspects of the problem of reconstructing objects from radiographs, Bulletin of the American Mathematical Society 83(6) (1977), 1227–1270.

Natterer

, The mathematics of computerized tomography, 2nd ed. Society for Industrial and Applied Mathematics, Philadelphia (1986).

Louis

, Incomplete data problems in X-ray computerized tomography 1. singular value decomposition of the limited angle transform, Numerische Mathematik 48 (1986), 251–262.

Frikel

and Quinto

, Characterization and reduction of artifacts in limited angle tomography, Inverse Problems 29(12) (2013), 1–21.

Quinto

, Artifacts and Visible Singularities in Limited Data X-Ray Tomography, Sensing & Imaging 18(1) (2017), 1–14.

Sidky

, Kao

, et al., Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, Journal of X-Ray Science and Technology 14(2) (2009), 119–139.

and Wang

, Compressed sensing based interior tomography, Physics in Medicine and Biology 54(9) (2009), 2791–2805.

Chen

, Jin

, et al., A limited-angle CT reconstruction method based on anisotropic tv minimization, Physics in Medicine and Biology 58(7) (2013), 2119–2141.

, Zhao

, et al., An image reconstruction model regularized by edge-preserving diffusion and smoothing for limited-angle computed tomography, Inverse Problems 35(8) (2019), 1–34.

10.

Hammernik

, Würfl

, et al., A deep learning architecture for limited-angle computed tomography reconstruction, Bildverarbeitung für die Medizin 12(14) (2017), 92–97.

11.

Huang

, Lu

, et al., Traditional machine learning techniques for streak artifact reduction in limited angle tomography, Bildverarbeitung für die Medizin 2018 11(13) (2018), 222–227.

12.

Bubba

, Kutyniok

, et al., Learning the invisible: a hybrid deep learning-shearlet framework for limited angle computed tomography, Inverse Problems 35(6) (2019), 1–38.

13.

Zhou

, Zhou

, et al., Limited view tomographic reconstruction using a cascaded residual dense spatial-channel attention network with projection data fidelity layer, IEEE Transactions on Medical Imaging 40(7) (2021), 1792–1804.

14.

Zhang

, Southern Medical University, China, Projection Data Optimization and Correction of Digital Breast Tomosynthesis Imaging (2017).

15.

, Cai

, et al., Promising generative adversarial network based sinogram inpainting method for ultra-limited-angle computed tomography imaging, Sensors 19(18) (2019), 1–21.

16.

Zhang

, Chen

, et al., Directional-TV algorithm for image reconstruction from limited-angular-range data. Medical Image Analysis 70 (2021), 1–15.

17.

Shi

and Qu

, A preconditioned Landweber iteration scheme for the limited-angle image reconstruction, Journal of X-Ray Science and Technology 29(6) (2021), 1045–1063.

18.

Wilkinson

, Rounding errors in algebraic processes, Prentice-Hall, Inc., Englewood Cliffs, USA (1963).

19.

Zhou

, Bovik

, Sheikh

and Simoncelli

, Image Qualifty Assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing 13(4) (2004), 600–612.