l 0 regularization based on a prior image incorporated non-local means for limited-angle X-ray CT reconstruction

Abstract

PURPOSES:

Restricted by the scanning environment in some CT imaging modalities, the acquired projection data are usually incomplete, which may lead to a limited-angle reconstruction problem. Thus, image quality usually suffers from the slope artifacts. The objective of this study is to first investigate the distorted domains of the reconstructed images which encounter the slope artifacts and then present a new iterative reconstruction method to address the limited-angle X-ray CT reconstruction problem.

METHODS:

The presented framework of new method exploits the structural similarity between the prior image and the reconstructed image aiming to compensate the distorted edges. Specifically, the new method utilizes l₀ regularization and wavelet tight framelets to suppress the slope artifacts and pursue the sparsity. New method includes following 4 steps to (1) address the data fidelity using SART; (2) compensate for the slope artifacts due to the missed projection data using the prior image and modified nonlocal means (PNLM); (3) utilize l₀ regularization to suppress the slope artifacts and pursue the sparsity of wavelet coefficients of the transformed image by using iterative hard thresholding (l₀W); and (4) apply an inverse wavelet transform to reconstruct image. In summary, this method is referred to as “l₀W-PNLM”.

RESULTS:

Numerical implementations showed that the presented l₀W-PNLM was superior to suppress the slope artifacts while preserving the edges of some features as compared to the commercial and other popular investigative algorithms. When the image to be reconstructed is inconsistent with the prior image, the new method can avoid or minimize the distorted edges in the reconstructed images. Quantitative assessments also showed that applying the new method obtained the highest image quality comparing to the existing algorithms.

CONCLUSIONS:

This study demonstrated that the presented l₀W-PNLM yielded higher image quality due to a number of unique characteristics, which include that (1) it utilizes the structural similarity between the reconstructed image and prior image to modify the distorted edges by slope artifacts; (2) it adopts wavelet tight frames to obtain the first and high derivative in several directions and levels; and (3) it takes advantage of l₀ regularization to promote the sparsity of wavelet coefficients, which is effective for the inhibition of the slope artifacts. Therefore, the new method can address the limited-angle CT reconstruction problem effectively and have practical significance.

Keywords

Computed tomography image reconstruction nonlocal means wavelet tight framelet prior image

1 Introduction

X-ray computed tomography (CT) is an effective modality to estimate the attenuation distribution of the scanned object, which is extensively applied in clinical diagnosis and non-destructive testing [1]. Generally, the obtained projection data are complete in classical CT application, for instance, the scanning angle needs a consecutive 180° for parallel reconstruction, or at least 180°+fan-angle for fan-beam reconstruction [2]. For specific cases, complete projection data are not efficient to use or even impossible to obtain, such as restricted scanning [3], dental CT [4], C-arm CT [5 –7], etc. In these situations, the obtained projection data may be limited-angle. Reconstructing from limited-angle projection data is more challenging than complete projection data. That’s because it might not be possible to accurately reconstruct the features and details of the reconstructed image from limited-angle projection data [8]. Meanwhile, the conventional analytic algorithm filtered backprojection (FBP) [9] or the simultaneous algebraic reconstruction technique (SART) [10] will not work for dealing with the limited-angle reconstruction [11], which leads to heavy distorted edges (slope artifacts) shown as following Fig. 1.

Fig.1

(a) Reference image and limited-angle scanning geometry in angular scope [180°, 300°], the red line of the object represents the point where the rays are not tangent to the edges. (b) Reconstructed image using FBP corresponding to the scanning trajectory in (a). (c) The difference image between reference image and reconstructed image.

Although it is difficult to settle limited-angle reconstruction problem, there exists analysis method for slope artifacts which are tangent to the boundary of some features in the object [8]. In reference [12], the slope artifacts which possess some directional property and are globally distributed have been researched. To surmount the slope artifacts, several efforts have been made to solve these thorny problems. Motivated by the success of the total variation (TV) algorithm applying in few-view reconstruction [13], the projection onto convex sets (POCS) with the total variation (TV) minimization as a representative was developed for limited-angle reconstruction [14]. Although TV-based algorithms have the ability to handle the localized artifacts and disturbances, they are not successful in revising the globalized artifacts due to the incomplete limited-angle projection data. Recently, l₀ regularization of wavelet coefficients based iterative reconstruction was developed to preferably suppress the slope artifacts and preserve the edges [15], but it will not work very well when the reconstructed image contains a large number of features and details.

The previous regularization methods are based on the prior condition of the piecewise constant of the reconstructed image. Subsequently, incorporating a previous image (prior image, for example, the reconstructed image from complete projection data or full-dose projection data) into the current image reconstruction framework is a promising research [16 –19]. A prior image constrained compressed sensing (PICCS) algorithm [17] as a representative was proposed to settle few-view reconstruction. However, it requires that the prior and current images should be taken at the same geometrical coordinates and very consistent. Recently, one image reconstruction algorithm [19] based on wavelet tight frame and prior image was proposed for limited-angle projection data, which is named as ‘l₀W-PI’ in this paper. It needs that the reconstructed image is consistent to the prior image, so it has a similar shortcoming with PICCS. With an edge-preserving nonlocal (NL) prior which takes the basic structure property of the image, nonlocal regularization for limited-angle X-ray CT was proposed [20, 21]. This idea was originally derived from image denoising [22], which considers the self-similarity properties in denoised image. Nonlocal means (NLM) filtering has a significant performance in image denoising with preserving the structural edges of reconstructed image. Recently, an iterative reconstruction using prior image induced nonlocal regularization was proposed for low-dose X-ray CT [23], where it takes a high-quality normal-dose scanning as a prior image.

The way to mitigate the slope artifacts for limited-angle CT, one side is to suppress the noise, the other side is to preserve the edges and details. Depending on limited-angle reconstruction problem, it might not be possible to precisely reconstruct the features of the image to be reconstructed [8]. Introducing the prior image into the limited-angle problem is an effective way. Whereas, in most cases, the prior image is not exactly consistent with the image to be reconstructed. Nevertheless, they are similar in structure.

For the sake of addressing the limited-angle reconstruction problem, motivated by the merits of l₀ regularization, wavelet tight framelets and NLM, in this work, we presented an iterative reconstruction method based on the structural similarity between prior and reconstructed images induced and modified NLM, and it was addressed via l₀ regularization of wavelet coefficients of the reconstructed image, which can be referred to as “l₀W-PNLM”. Specifically, on one hand, the presented framework exploits the structural similarity between prior image and reconstructed image, which can compensate the distorted edges; on the other hand, this work utilizes l₀ regularization and wavelet tight framelets aiming to suppress the slope artifacts.

The rest of this work is organized as follows: in Section II, compared algorithms are prepared and the method in this paper is presented; in Section III, numerical experiments are implemented to verify the feasibility and effectiveness of the presented l₀W-PNLM; in Section IV, conclusion and discussion are given.

2 Methods

2.1 Review of TV, PICCS and l₀W-PI algorithms

Total variation (TV) and prior image constrained compressed sensing (PICCS) are both essentially iterative algorithms to reconstruct an unknown image x via minimizing l₁ regularization based optimization problems. In conventional TV algorithm, there is not a prior image used directly, while on assumption that the desired image under some particular transformation is sparse. From a mathematical point of view, let x denote the image to be reconstructed, and TV regularization based optimization problem [13] can be expressed as following: $min_{x} {∥ x ∥}_{TV}, s . t . Ax = b, x \geq 0 .$ (1)

Where A is a system matrix, b denotes the obtained projection data, and ${∥ x ∥}_{TV} = \sum_{i, j} \sqrt{{(x_{i, j} - x_{i - 1, j})}^{2} + {(x_{i, j} - x_{i, j - 1})}^{2}}$ . TV algorithm is efficiently performed with incomplete projection data for fan-beam CT, especially for few-view situation. Whereas, TV algorithm will not perform well for the situation that is reconstructing the image which contains lots of the features and details (such as vivo animal images) under seriously under-sampled [24, 25]. Meanwhile, it is difficult to manipulate the slope artifacts due to lack of enough projection data.

Aiming to enable the accurate reconstruction from the seriously under-sampled projection data in dynamic CT, a new framework PICCS was proposed [17, 26], whose idea is to take advantage of the relationship between the prior image and the reconstructed image. Let x_p denote the prior image, and then the associated optimization problem can be represented as following: $min_{x} {μ {∥ x - x_{p} ∥}_{TV} + (1 - μ) {∥ x ∥}_{TV}}, s . t . Ax = b, x \geq 0 .$ (2)

Where the balancing parameter μ is set to [0,1]. Supposing that μ = 0, then the problem (2) is equivalent to the problem (1). In the PICCS objective function, on one hand, it considers the sparsity of the gradient transformation of the image x itself, on the other hand, it utilizes the similarity and consistency between the prior image x_p and the reconstructed image x. This enables accurate image reconstruction from few-view projection data. For the limited-angle CT reconstruction problem, although PICCS algorithm will work well when the prior image is similar and consistent to the reconstructed image, it will fail when the difference between prior and reconstructed images is not very small.

Considering that the reconstructed image can be fine approximated using sparse wavelet coefficients under some wavelet tight frames, to deal with limited-angle CT reconstruction problem, a minimization model based on wavelet tight frames and a prior image was proposed [19] shown as follows. It is denoted as ‘l₀W-PI’ in this paper. $min_{x} {\frac{1}{2} {∥ Ax - b ∥}_{2}^{2} + \frac{γ}{2} {∥ {(Wx)}_{G} - {({Wx}_{p})}_{G} ∥}_{2}^{2} + β {∥ {(Wx)}_{L} ∥}_{0}}, s . t . x \geq 0 .$ (3)

Where W shows the wavelet tight frames; Wx = ((Wx) _L, (Wx) _G), (Wx) _L denotes the low-frequency of the image x under W transform, (Wx) _G and (Wx_p) _G denote the high-frequency of the images x and x_p under W transform, respectively; ∥ · ∥ ₀ represents the number of nonzero terms; γ and β are the positive regularization parameters. l₀W-PI algorithm can obtain good quality image when the reconstructed image is similar to the prior image, but it exists the same shortcoming with PICCS, that is, it will lead to the distorted edges in the reconstructed results where the reconstructed image is different from the prior image. The detailed solution process can be referred to the reference [19].

2.2 Modified NLM

Nonlocal means (NLM) was originally used for image denoising [22]. It amounts to mean that the denoised image x at point n_i (Denoted as NLM (x) (n_i)) is a weighted mean of the values for all points which belongs to its search-window Ω_x. x (n_i) denotes the intensity of the image x at the point n_i. Its simple formula can be defined as following: $NLM (x) (n_{i}) = \frac{1}{c (x (n_{i}))} \sum_{n_{j} \in Ω_{x}} w (n_{i}, n_{j}) x (n_{j}) .$ (4)

Where V_i, V_j ⊂ Ω_x are the position matrixes of the image which denote the two patch-windows with the same size, and their center points are n_i and n_j, respectively. x (V_i) ⊂ x and x (V_j) ⊂ x denote the intensity of the image x at the patch-windows V_i and V_j respectively. The correlation coefficient $w (n_{i}, n_{j}) = exp {- \frac{{∥ x (V_{i}) - x (V_{j}) ∥}_{2}^{2}}{h^{2}}}$ , and normalized factor c (x (n_i)) = ∑_{n_j∈Ω_x}w (n_i, n_j). h is a smoothing parameter which is usually proportional to standard deviation of the image noise. In general, the situation is that it handles the image x itself like the reference [22] or [27], which considers the local similarity of the images.

In this work, we modified the nonlocal means filtering to promote prior image based limited-angle reconstruction. Let the point n_i belong to reconstructed image x, and its search-window is Ω_x and its patch-window is V_i. Let the point n_j belong to prior image x_p, and its relevant symbols are Ω_{x
_p} and V_j. Then, the modified nonlocal means incorporated with the prior image (PNLM) of reconstructed image x at the point n_i can be expressed as: $PNLM (x) (n_{i}) = \frac{1}{c_{p} (x_{p} (n_{j}))} \sum_{n_{j} \in Ω_{x_{p}}} w_{p} (n_{i}, n_{j}) x_{p} (n_{j}) .$ (5)

Different from Equation (4), the correlation coefficient $w_{p} (n_{i}, n_{j}) = exp {- \frac{{∥ x (V_{i}) - x_{p} (V_{j}) ∥}_{2}^{2}}{h^{2}}}$ and normalized factor c_p (x_p (n_j)) = ∑_{n_j∈Ω_{x
_p}}w_p (n_i, n_j); n_i ∈ Ω_x belongs to reconstructed image x and n_j ∈ Ω_{x
_p} belongs to prior image x_p in the Equation (5). This considers the structural similarity between reconstructed image and prior image. Since the slope artifacts are distributed globally in the reconstructed results like Fig. 1, the main impact is associated with missed angles. Hence, Equation (5) is utilized to amend the edges of the reconstructed image in heavy slope artifact location as shown in Fig. 2.

Fig.2

It shows the modified NLM. (a) represents the reconstructed result using FBP. (b) shows the prior image.

2.3 Wavelet tight frames

When the obtained projection data are seriously limited-angle, it will lead to the slope artifacts in reconstructed image. In order to surmount the slope artifacts for limit-angle CT reconstruction, sparse transforms, such as gradient, Haar transform [28], S-transform [29], and wavelet frames [30 –33], are necessary to be adopted to handle the reconstructed image. In consideration of several advantages of the piecewise-constant linear B-spline framelet transform as follows: on one side, it can sparsely represent the piecewise smooth functions, such as images [30]; on the other side, it can be seen as a generalization of total variation [32]; additionally, it can be treated as an extension of the wavelet which provides more multilevel redundancy in sparsely representing the image [32]. In this work, the piece-wise constant B-spline frames W (its inverse transform is W^T) was utilized which can be constructed by the reference [30]. Its relevant filters can be represented as: $\frac{1}{4} [1, 2, 1], \frac{\sqrt{2}}{4} [1, 0, - 1], \frac{1}{4} [- 1, 2, - 1] .$

Its detailed theory can refer to the references [30, 32].

The wavelet tight frames are utilized to transform the image which contains slope artifacts. When the wavelet coefficients α are obtained, the following model [34] will be adopted to deal with that in this paper. $min_{α} {{∥ α - \tilde{α} ∥}_{2}^{2} + λ {∥ α ∥}_{0}} .$ (6)

The model (6) considers all the coefficients. The reason is that the slope artifacts are not only in the high frequency components of the wavelet decomposition, but also in low frequency. In addition, l₀ regularization has the capability of addressing the slope artifacts.

2.4 Presented reconstruction method

In consideration of the existing shortcomings of TV, PICCS and l₀W-PI algorithms for settling the limited-angle reconstruction, we presented an iterative reconstruction method based on the prior image induced and modified NLM and l₀ regularization of wavelet coefficients of the reconstructed image, which is referred to as “l₀W-PNLM”. Subsequently, the detailed process of the presented method is introduced. It includes four steps as following:

In the first step, SART technique [10] is utilized to reconstruct the image from the limited-angle projection data with noise. Rearrange the image x to be a vector, and it can be denoted as: $y_{j}^{n} = x_{j}^{n} + ω \frac{1}{\sum_{i = 1}^{M} a_{ij}} \sum_{i = 1}^{M} \frac{a_{ij}}{\sum_{j = 1}^{N} a_{ij}} (b_{i} - A_{i} x^{n}), j = 1, 2, \dots, N .$ (7)

Where ω is a step length, A ∈ R^MN is a system matrix, A_i denotes the ith row of system matrix A, a_ij is the ith row and jth column element of A, b_i represents the ith component of projection data b, M is the total number of X-rays and N denotes the total number of image pixels. In brief, it can be denoted as yⁿ = SART (xⁿ).

In the second step, when iteration number equals to a certain number, then Equation (5) is used to modify the edges by the damage of slope artifacts. Rearrange the result yⁿ to be a matrix. Let Ω be the whole domain of the reconstructed image and Ω_x be the search-window. Ω_x ⊂ Ω and its formation can be expressed as following: ${\tilde{y}}^{n} (n_{i}) = {\begin{matrix} PNLM (y^{n}) (n_{i}), n_{i} \in Ω_{x}; \\ y^{n} (n_{i}), n_{i} \notin Ω_{x} . \end{matrix} n_{i} \in Ω .$ (8)

In the third step, the wavelet tight frames are adopted to transform ${\tilde{y}}^{n}$ , denoted by ${\tilde{α}}^{n} = W ({\tilde{y}}^{n})$ . Let λ be a parameter and $Ω_{W}^{kl}$ a wavelet domain which represents the kth of the lth layer. L denotes the total number of the wavelet transform layers and K is the total number of wavelet coefficients of each layer. Then, to solve the model (6), hard thresholding algorithm [34] is exploited to deal with the wavelet coefficients ${\tilde{α}}^{n}$ , denoted by $α^{n + 1} (n_{i}) = {\begin{matrix} {\tilde{α}}^{n} (n_{i}), | {\tilde{α}}^{n} (n_{i}) | > \sqrt{λ}; \\ 0, | {\tilde{α}}^{n} (n_{i}) | \leq \sqrt{λ} . \end{matrix} n_{i} \in Ω_{W}^{kl}, k = 1, \dots, K, l = 1, \dots, L .$ (9)

In the forth step, do the inverse wavelet tight frames transformation according to αⁿ⁺¹, which can be expressed as xⁿ⁺¹ = W^T (αⁿ⁺¹).

Summarize the above all steps and let N_max be the maximum iteration number. Then the pseudo-code of the presented l₀W-PNLM can be described as shown in Table 1. One algorithm only includes Step 1 and Step 2, which is named as ‘PNLM’, and the other algorithm only includes Step 1, Step 3 and Step 4, which is named as ‘l₀W’ for the following comparative experiments.

Table 1

The pseudo-code of the presented l₀W-PNLM

Initialization: n = 1, ω = 1, K = 9, L = 4, x₁ = 0, K α¹ = Wx¹, h, λ, V_i, V_j, Ω_x, Ω_{x
_p}, N_i1, N_i2, N_max.

While (n ≤ N_max) do

Step 1: If n < N_i1 do y_n = SART (x_n), else y_n = TV (x_n) end (if);

Step 2: If n = N_i2 do

{\tilde{y}}^{n} (n_{i}) = {\begin{matrix} PNLM (y^{n}) (n_{i}), n_{i} \in Ω_{x}; \\ y^{n} (n_{i}), n_{i} \notin Ω_{x} . \end{matrix} n_{i} \in Ω

{\tilde{y}}^{n} : = 0 . 3^{*} y^{n} + 0 . 7^{*} {\tilde{y}}^{n} .

end (if);

Step 3:

{\tilde{α}}^{n} = W ({\tilde{y}}^{n})

α^{n + 1} (n_{i}) = {\begin{matrix} {\tilde{α}}^{n} (n_{i}), | {\tilde{α}}^{n} (n_{i}) | > \sqrt{λ}; \\ 0, | {\tilde{α}}^{n} (n_{i}) | \leq \sqrt{λ} . \end{matrix} n_{i} \in Ω_{W}^{kl}, k = 1, \dots, K, l = 1, \dots, L .

Step 4: xⁿ⁺¹ = W^T (αⁿ⁺¹), n:=n+1.

end (While)

3 Numerical implementations

In this section, numerical experiments for limited-angle reconstruction are implemented to validate and evaluate the effectiveness of the presented l₀W-PNLM, compared with TV and PICCS algorithms. In the experiments, quantitative assessments are investigated to verify the performance of l₀W-PNLM. The evaluation indexes include the root mean squared error (RMSE) [35], peak signal to noise ratio (PSNR) [35], universal image quality index (UQI) [36] and structural similarity index (SSIM) [37] shown as following: $RMSE = \sqrt{\frac{1}{N} \sum_{n_{i}} {(x (n_{i}) - x^{r} (n_{i}))}^{2}},$ (10) $PSNR = 10 {log}_{10} \frac{{∥ x ∥}_{\infty}^{2}}{{RMSE}^{2}} (dB) .$ (11) $UQI = \frac{4 μ_{x} μ_{x^{r}} Cov {x, x^{r}}}{(σ_{x}^{2} + σ_{x^{r}}^{2}) (μ_{x}^{2} + μ_{x^{r}}^{2})} .$ (12) $SSIM = \frac{2 μ_{x} μ_{x^{r}} (2 Cov {x, x^{r}} + c_{2})}{(μ_{x}^{2} + μ_{x^{r}}^{2} + c_{1}) (σ_{x}^{2} + σ_{x^{r}}^{2} + c_{2})} .$ (13)

Where x (n_i) denotes the intensity at the point n_i of reconstructed image x; x^r is the reference image and x^r (n_i) is the intensity at the point n_i of x^r; N is the total number of the image pixels; ${∥ x^{r} ∥}_{\infty} = max_{n_{i}} | x^{r} (n_{i}) |$ ; μ_x and μ_{x
^r} are the mean values of x and x^r; σ_x and σ_{x
^r} are the standard deviations of x and x^r; Cov{ x, x^r } is the co-variance between x and x^r; c₁ = (0.01R) ² and c₂ = (0.03R) ², R is the dynamic range of pixel values.

In the experiments, the parameters of our method are chose by trial and error for better reconstructed image quality. All the experiments are implemented on 2.40 GHz intel(R) Core(TM) i5-4210U CPU processor with 4G memory and coded using Microsoft Visual Studio 2010 C++ and Matlab R2013a.

3.1 Simulated NCAT experiments

In the experiments, a digital NURBS based cardiac-torso (NCAT) phantom [38] as shown in Fig. 3(a) was tested, where its four corners are the local zoom-in figures corresponding to the marked numbers, respectively. Supposing that prior image and reference image are registered, then in Fig. 3(b) represents the prior image and (c) shows the difference image between (a) and (b). From Fig. 3(c), it can be observed that the structure of prior image is similar to reference image, while the difference is not small, such as the place where the central part of the prior image is smaller than that of the reference image. It follows that some defects or details of the image to be reconstructed will not be covered by the prior image. The simulated projection data were generated via projecting a 256×256 discretized NCAT phantom with adding Gaussian noise whose mean value is zero and standard deviation is 0.1% ∥b ∥ _∞. The geometry scanning parameters for simulated CT imaging system are listed in Table 2. Three scanning angular scopes are [0, 90°], [0, 100°] and [0, 120°].

Fig.3

It shows the reference NCAT image, prior image and difference image between two former. Four corners of (a) are the local zoom-in figures corresponding to the marked numbers, respectively. The gray scale window is [0,255].

Table 2

Geometry scanning parameters of NCAT for simulated CT imaging system

Parameter	Value
The distance from X-ray source to rotation center	900.0mm
The distance from X-ray source to detector	1300.0mm
The diameter of field of view	255.0mm
Interval angle between two adjacent projection views	0.6679°
The number of detector units	256
Image size	256×256
Pixel size of the image	1.0×1.0mm²

For the experiments, the parameters are determined by trial and error for better reconstructed image quality. For all algorithms, the maximal iteration number N_max is set to 800. Another parameters are determined as follows:

In TV algorithm, two parameters are chosen: one is the parameter of SART step, which is set to 1; the other is in TV step, which is set to 0.2. In PICCS algorithm, the parameter of fidelity term is set to 10 and μ is set to 0.2. Similarly in SART step, the step parameter is set to 1.

For the presented l₀W-PNLM, N_i1 = 550 and N_i2 < 600 is set to an iteration number which can be divisible by 100. In Step 1, ω = 1 and the parameter for TV step is set to 0.06. In Step 2, the size of search-windows Ω_x and Ω_{x
_p} should be enough to obtain more redundant information, but the large size will increase the computational burden of the algorithm. In this paper, taking the prior knowledge that the limited-angle artifacts are affected by the missing projection data into account, search-windows Ω_x and Ω_{x
_p} are both set to [7 128]×[128 240] which are the coordinates of the image x, [7 128] represents the rows 7–128 of the image x, and [128 240] denotes the columns 128–240 of the image x. The size of patch-windows V_i and V_j is set to 5×5. For smoothing parameter h, previous studies show that it is proportional to standard deviation of the image noise, and it is set to 0.05 in this paper. In Step 3, the parameter of hard thresholding λ is determined to 0.0035² via balancing data consistency and sparsity.

3.2 NCAT results and quantitative assessments

Figure 4 shows the results reconstructed using TV, PICCS and l₀W-PNLM from incomplete projection data in angular scopes [0, 120°], [0, 100°] and [0, 90°], respectively. In order to better show the contrast details, from Figs. 5 to 8, they are the local zoom-in images corresponding to the marked number 1 to 4. It can be observed that the edges (especially in Figs. 6, 8, 9) of images reconstructed using TV algorithm are distorted due to the missing angular projection data like Fig. 1. It can be still seen that the results using PICCS algorithm has distorted edges especially in some places where reconstructed image is different from prior image. Compared with TV and PICCS algorithms, the presented l₀W-PNLM has more ability to suppress the slope artifacts and preserve the edges. Although the scanning angular scope decreases, our method can still maintain the good image quality.

Fig.4

From left column to right column, it shows the results reconstructed using TV, PICCS and l₀W-PNLM algorithms from limited-angle projection data in angular scopes [0,120°], [0,100°] and [0,90°], respectively. The gray scale window is [0,255].

Fig.5

It presents the local zoom-in images of Fig. 4 corresponding to the marked number 1. The gray scale window is [0,255].

Fig.6

It presents the local zoom-in images of Fig. 4 corresponding to the marked number 2. The gray scale window is [0,255].

Fig.7

It presents the local zoom-in images of Fig. 4 corresponding to the marked number 3. The gray scale window is [0,255].

Fig.8

It presents the local zoom-in images of Fig. 4 corresponding to the marked number 4. The gray scale window is [0,255].

Fig.9

It presents the quantitative assessments of the images reconstructed using TV, PICCS and l₀W-PNLM. (a) shows the RMSE. (b) shows the PSNR. (c) shows the UQI. (d) shows the SSIM.

In addition to visual comparison, quantitative assessments including RMSE, PSNR, UQI and SSIM are implemented to show the differences among TV, PICCS and l₀W-PNLM algorithms. Figure 9 shows the assessments of the images reconstructed using TV, PICCS and l₀W-PNLM algorithms. From Fig. 9(a), it can be seen that our method l₀W-PNLM can get better indicators, which is as the angular scope is decreasing, l₀W-PNLM can maintain the smallest RMSE value among the compared three algorithms. From Fig. 9(b), l₀W-PNLM can obtain the largest PSNR value among them. Figure 9(c) and (d) verify the same points. In these results, PICCS algorithm is not very effective for the issue that the reconstructed image and prior image are not very similar. Since TV algorithm does not utilize prior image, its effect is not obvious for the severely limited-angle reconstruction.

3.3 Simulated Chest experiments

To further perform the effectiveness of the presented method, a Chest phantom [39] as shown in Fig. 10(a) is tested. Figure 10(b) is the prior image and (c) represents the difference image between the reference image and prior image. Figure 10(c) shows that the structure of (a) is similar to (c), but not consistent. The simulated projection data were generated via projecting a 256×256 discretized Chest phantom with adding Gaussian noise whose mean value is zero and standard deviation is 0.1% ∥b ∥ _∞. Three scanning angular scopes are [0,90°], [0,100°] and [0,120°]. The simulated geometrical scanning parameters for limited-angle CT are the same with NCAT experiments listed in Table 2 except the pixel size which is set to 0.625×0.625mm².

Fig.10

It shows the reference Chest image, prior image and difference image between two former. Above two corners of (a) are the local zoom-in images corresponding to the marked red rectangles, respectively. The gray scale window is [0,255].

For all algorithms, the maximal iteration number N_max is set to 800 except 200 for PNLM algorithm. Another parameters are determined as follows:

In TV algorithm, two parameters are chosen: one is the parameter of SART step, which is set to 1; the other is in TV step, which is set to 0.16. In PICCS algorithm, the parameter of fidelity term is set to 100 and μ is set to 0.2. Similarly in SART step, the step parameter is set to 1. In l₀W-PI algorithm, the main parameters γ = 30 and β = 0.002, and the other parameters can be referred to the reference [19].

For the presented l₀W-PNLM, N_i1 = 680 and N_i2 < 200 is set to an iteration number which can be divisible by 60. In Step 1, ω is set to 1, 0.6 and 0.5 for the angular scope [0,120°], [0,100°] and [0, 90°], respectively. And the parameter for TV step is all set to 0.06. In Step 2, the search-windows Ω_x and Ω_{x
_p} are both set to [7 140]×[109 240] which are the coordinates of the image x, [7 140] represents the rows 7–140 of the image x, and [109 240] denotes the columns 109–240 of the image x. The size of patch-windows V_i and V_j is set to 5×5. For smoothing parameter h, previous studies show that it is proportional to standard deviation of the image noise, and it is set to 0.05 in this paper. In Step 3, the parameter of hard thresholding λ is determined to 0.0047² via balancing data consistency and sparsity.

For PNLM algorithm, in Step 1, ω = 1 for SART step. The parameter for TV step is set to 0.2 for angular scope [0, 120°], and 0.06 for angular scopes [0, 100°] and [0, 90°]. The other corresponding parameters are the same with l₀W-PNLM.

3.4 Chest results and quantitative assessments

Figures 11 and 12 show the results reconstructed using TV, PICCS, l₀W-PI, PNLM, l₀W and l₀W-PNLM from incomplete projection data in angular scopes [0, 120°], [0, 100°] and [0, 90°], respectively. To better exhibit the details, the corresponding zoom-in images are placed in the above two corners in each reconstructed image. It can be observed that although PICCS and l₀W-PI algorithms both have the ability to deal with slope artifacts, some edges are distorted as shown in the red arrow marks. They are not very effective for the issue that the reconstructed image and prior image are not very similar. PNLM has some effect in dealing with distorted edges, since it utilizes the structural similarity between the reconstructed image and prior image, but not work well for the heavy slope artifacts. TV algorithm leads to the distorted edges in reconstructed results, while l₀W algorithm can better control the slope artifacts. By contrast, l₀W-PNLM has more capability of addressing the limited-angle reconstruction problem. As the angular scope decreasing, l₀W-PNLM can maintain good image quality and preserve the edges.

Fig.11

It shows the Chest results. From the first row to the third row, it represents the results reconstructed using TV, PICCS and l₀W-PI algorithms from limited-angle projection data in angular scopes [0, 120°], [0, 100°] and [0, 90°], respectively. Above two corners of each image are the zoom-in images of the marked red rectangles in Fig. 10 (a). The gray scale window is [0,255].

Fig.12

It shows the Chest results. From the first row to the third row, it represents the results reconstructed using PNLM, l₀W and l₀W-PNLM algorithms from limited-angle projection data in angular scopes [0, 120°], [0, 100°] and [0, 90°], respectively. Above two corners of each image are the zoom-in images of the marked red rectangles in Fig. 10 (a). The gray scale window is [0,255].

Figure 13 represents the quantitative assessments of the Chest results. It can be seen that l₀W-PNLM obtains the best indicators and PNLM takes the second place. Although l₀W algorithm cannot obtain the overall performance, it gets the high SSIM indicator. It indicates that l₀W algorithm has the ability to reconstruct the structure of the image. PICCS and l₀W-PI have the weak indicators. TV gets the worse one.

Fig.13

It presents the quantitative assessments of the images reconstructed using TV, PICCS, l₀W-PI, PNLM, l₀W and l₀W-PNLM. (a) shows the RMSE. (b) shows the PSNR. (c) shows the UQI. (d) shows the SSIM.

Compared to TV, PICCS, l₀W-PI, PNLM and l₀W algorithms, what are the advantages of the presented l₀W-PNLM? On one aspect, it utilizes the structural similarity between the reconstructed image and prior image, which can modify the distorted edges suffering from slope artifacts. Additionally, this method adopts wavelet tight frames which can obtain the first and high derivative in several directions and levels. Moreover, it takes advantage of l₀ regularization to promote the sparsity of wavelet coefficients, which is effective for the inhibition of the slope artifacts. And furthermore, the incorporated prior image does not cover up some details and features of the reconstructed image. The presented method carries forward the advantages of the above algorithms.

4 Discussion

In this article, we presented an iterative reconstruction method l₀W-PNLM which is based on the structural similarity between prior image and reconstructed image induced modified NLM, and it is addressed via l₀ regularization of wavelet coefficients of the reconstructed image for limited-angle CT reconstruction. It not only takes advantage of the structural similarity between reconstructed image and prior image via modified NLM, but also adopts the advantages of suppressing the slope artifacts via l₀ regularization and wavelet tight framelets. We have evaluated and demonstrated the performance of l₀W-PNLM in addressing limited-angle reconstruction problem. As numerical implementations indicate that the presented method can yield better reconstructed images and quantitative assessments compared to TV, PICCS, l₀W-PI, PNLM and l₀W algorithms.

There might be some aspects of the presented method that are relevant and useful for medical and industrial limited-angle CT imaging. One aspect, the sparse hypothesis of the image under some transformation is quite reasonable for many objects in medical and industrial applications. The other aspect, for medical imaging, the prior image can be obtained in the patient’s previous examination and for industrial imaging, the prior image can be acquired before using. Therefore, it has some practical significance. However, the parameter selection of the presented method is by trial and error and is not adaptive. In the future, we will research the adaptivity of these parameters, and improve the presented method further for addressing limited-angle X-ray CT reconstruction problem.

Footnotes

Acknowledgment

This work is supported by the National Natural Science Foundation of China (61771003) and National Instrumentation Program of China (2013YQ030629).

References

Mihailidis

D.P.D.

, Computed tomography from photon statistics to modern cone-beam CT, Medical Physics 36(8) (2009), 3858–3858.

Tam

K.C.

, Reducing the fan-beam scanning angular range, Physics in Medicine and Biology 33(8) (1988), 955–967.

Gao

, Zhang

, Xing

et al., Volumetric imaging from a multisegment straight-line trajectory and a practical reconstruction algorithm, Optical Engineering 46(7) (2007), 077004.

Cho

M.K.

, Kim

H.K.

, Youn

et al., A feasibility study of digital tomosynthesis for volumetric dental imaging, Journal of Instrumentation 7(3) (2012), 55–61.

, Noo

, Lauritsch

et al., Extended Volume Image Reconstruction Using the Ellipse-Line-Ellipse Trajectory for a C-arm System, In Proceedings of the 12th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine (2013), 245–248.

Kuo

J.C.

, Lu

N.H.

, Huang

Y.H.

et al., Analyzed outcome for embolization of hepatic carcinoma by C-Arm angiographic computed tomography, Journal of X-ray Science and Technology 22(5) (2014), 645–651.

, Maier

, Lauritsch

et al., Axially extended-volume C-arm CT using a reverse helical trajectory in the interventional room, IEEE Transactions on Medical Imaging 34(1) (2015), 203–215.

Quinto

E.T.

, Artifacts and visible singularities in limited data X-Ray tomography, Sensing and Imaging 18(1) (2017), 9.

Kak

A.C.

and Slaney

, Principles of computerized tomographic imaging, Medical Physics 29(1) (2002), 107–107.

10.

Andersen

A.H.

and Kak

A.C.

, Simultaneous algebraic reconstruction technique (SART): A superior implementation of the art algorithm, Ultrasonic Imaging 6(1) (1984), 81.

11.

Natterer, Frank, Wubbeling, et al., Mathematical methods in image reconstruction, Medical Physics 29(1) (2001), 107–108.

12.

and Ye

J.C.

, Multi-scale wavelet domain residual learning for limited-angle CT reconstruction, (2017), arXiv:1703.01382.

13.

Sidky

E.Y.

, Kao

C.M.

and Pan

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

14.

Sidky

E.Y.

and Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in Medicine and Biology 53(17) (2008), 4777–4807.

15.

Zeng

and Wang

, Error bounds and stability in the l₀ regularized for CT reconstruction from small projections, Inverse Problems and Imaging 10(3) (2017), 829–853.

16.

Jianhua

Z.H.

, Gao

et al., Iterative image reconstruction for cerebral perfusion CT using a pre-contrast scan induced edge-preserving prior, Physics in Medicine & Biology 57(22) (2012), 7519.

17.

Chen

G.H.

, Theriault-Lauzier

, Tang

et al., Time-resolved interventional cardiac C-arm cone-beam CT: An application of the PICCS algorithm, IEEE Transactions on Medical Imaging 31(4) (2012), 907–923.

18.

Ouyang

, Solberg

and Wang

, Noise reduction in low-dose cone beam CT by incorporating prior volumetric image information, Medical Physics 39(5) (2012), 2569–2577.

19.

Wang

, Zeng

, Guo

et al., Wavelet tight frame and prior image-based image reconstruction from limited-angle projection data, Inverse Problems & Imaging 11(6) (2017), 917–948.

20.

D.F.

and Fessler

J.A.

, Edge-preserving tomographic reconstruction with nonlocal regularization, IEEE Transactions on Medical Imaging 21(2) (2002), 159–173.

21.

Liu

, Katsevich

and Yu

, Interior tomography with curvelet-based regularization, Journal of X-ray Science and Technology 25(1) (2017), 1–13.

22.

Buades

, Coll

and Morel

J.M.

, A non-local algorithm for image denoising, IEEE Computer Society Conference on Computer Vision and Pattern Recognition IEEE Computer Society 2(7) (2005), 60–65.

23.

Zhang

, Huang

, Ma

et al., Iterative reconstruction for X-Ray computed tomography using prior-image induced nonlocal regularization, IEEE Trans Biomed Eng 61(9) (2014), 2367.

24.

Tang

, Nett

B.E.

and Chen

G.H.

, Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms, Physics in Medicine and Biology 54(19) (2009), 5781.

25.

Block

K.T.

, Uecker

and Frahm

, Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint, Magnetic Resonance in Medicine 57(6) (2007), 1086–1098.

26.

Chen

G.H.

, Tang

and Leng

, Prior image constrained compressed sensing (PICCS), Medical Physics 35(2) (2008), 660–663.

27.

, Chen

, Wu

et al., Iterative image reconstruction using modified non-local means filtering for limited-angle computed tomography, Physica Medica 32(9) (2016), 1041–1051.

28.

Garduño

, Herman

G.T.

and Davidi

, Reconstruction from a few projections by l₁-minimization of the haar transform, Inverse Problems 27(5) (2011), 55006–55018(13).

29.

Luo

, Liu

, Li

et al., Image reconstruction from sparse projections using S-transform, Journal of Mathematical Imaging and Vision 43(3) (2012), 227–239.

30.

Ron

and Shen

, Affine systems in L2 (Rd): The analysis of the analysis operator, Journal of Functional Analysis 148(2) (1997), 408–447.

31.

Gao

, Cai

J.F.

, Shen

et al., Robust principal component analysis-based four-dimensional computed tomography, Physics in Medicine and Biology 56(11) (2011), 3181–3198.

32.

Zhao

, Gao

, Ding

et al., Tight-frame based iterative image reconstruction for spectral breast CT, Medical Physics 40(3) (2013), 031905.

33.

Dong

, Shen

and Xie

, Image restoration: A general wavelet frame based model and its asymptotic analysis, (2016), arXiv:1602.05332.

34.

Blumensath

and Davies

M.E.

, Iterative hard thresholding for compressed sensing, Applied and Computational Harmonic Analysis 27(3) (2008), 265–274.

35.

Storath

, Weinmann

, Frikel

et al., Joint image reconstruction and segmentation using the potts model, Inverse Problems 31(2) (2014), 25003–25031(29).

36.

Wang

and Bovik

A.C.

, A universal image quality index, IEEE Signal Processing Letters 9(3) (2002), 81–84.

37.

Wang

, Bovik

A.C.

, Sheikh

H.R.

and Simoncelli

E.P.

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

38.

Xcat phantom website Available: http://www.hopkinsradiology.org/DMIP/Research/xcat.

39.

Chest and abdomen phantom website. Available: http://www9.informatik.uni-erlangen.de/External/vollib/.

l 0 regularization based on a prior image incorporated non-local means for limited-angle X-ray CT reconstruction

Abstract

PURPOSES:

METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1 Introduction

2.1 Review of TV, PICCS and l0W-PI algorithms

Footnotes

Acknowledgment

References

2.1 Review of TV, PICCS and l₀W-PI algorithms