Total variation combining nonlocal means filtration for image reconstruction in X-ray computed tomography

Abstract

BACKGROUND:

Image reconstruction for realistic medical images under incomplete observation is still one of the core tasks for computed tomography (CT). However, the stair-case artifacts of Total variation (TV) based ones have restricted the usage of the reconstructed images.

OBJECTIVE:

This work aims to propose and test an accurate and efficient algorithm to improve reconstruction quality under the idea of synergy between local and nonlocal regularizations.

METHODS:

The total variation combining the nonlocal means filtration is proposed and the alternating direction method of multipliers is utilized to develop an efficient algorithm. The first order approximation of linear expansion at intermediate point is applied to overcome the computation of the huge CT system matrix.

RESULTS:

The proposed method improves root mean squared error by 25.6% compared to the recent block-matching sparsity regularization (BMSR) on simulation dataset of 19 views. The structure similarities of image of the new method is higher than 0.95, while that of BMSR is about 0.92. Moreover, on real rabbit dataset of 20 views, the peak signal-to-noise ratio (PSNR) of the new method is 36.84, while using other methods PSNR are lower than 35.81.

CONCLUSIONS:

The proposed method shows advantages on noise suppression and detail preservations over the competing algorithms used in CT image reconstruction.

Keywords

Computed tomography Image reconstruction nonlocal means filtration total variation alternating direction method of multipliers first order approximation

1 Introduction

Image reconstruction is one of the core tasks of computed tomography (CT) which is a kind of typical inverse problems [1]. In the discrete-to-discrete model [2], the data acquisition can be approximately described as Ax = b + r, where vector b is the observation data and it is usually contaminated by noise r, and A is the sensing matrix or system matrix which is determined by the scanning geometry. The core task of reconstruction is to find x according to scanning geometry (A) and the noisy observation (b + r), which is equivalent to solve ${find x \in ℝ^{n} | Ax = b + r}$ . When facing with incomplete data, the reconstruction problem often needs efficient regularization functional [3]. By exploiting the sparsity of the transformed domain of the intended image, the reconstruction can be modeled as a constrained optimization problem as x ∈ {min R (x) | Ax = b + r}, where R (x) is the regularization function applied on image x. The key problem to develop a useful algorithm of image reconstruction relies on a proper regularization model and an efficient solvation algorithm.

The image sparsity enabled the formation of regularization function [4]. With a given image, the gradient domain is highly sparse if it is piece-wise constant which driven the application of total variation (TV) utilized in the optimization. The algorithm of adaptive steepest descent-projection onto convex set [5] proposed by Sidky and Pan has drawn much attention on TV based methods, with its improved variants frequently proposed [6, 7]. Moreover, the alternating direction method of multipliers [8] (ADMM) introduced into CT image reconstruction on L1 norm based TV (TV-L1) [9] model by Zhang improved the efficiency of the convergence behavior. It must be pointed out that the TV-L1 is a kind of local model which only focus on neighbor pixels of image. When applying TV on medical images under incomplete observation, the so-called stair case blurring are often present as artifacts which severely deteriorate the image quality. Recent years, Kim has proposed the reweighted L1 norm based nonlocal TV (NLTV) [10], which is solved by simultaneous algebraic reconstruction technique [11] (SART) and a followed NLTV post-processing. Although NLTV following SART shows better compared to TV based ones, the stair-case artifacts are still not satisfied removed. More recently, the block-matching sparsity regularization [12, 13] (BMSR) on nonlocal pixels proposed by Cai shows advantages on more realistic medical images. Since the piece-wise constant property is often not applicable on medical images, BMSR shows advantages on suppression of stair case artifacts. In spite of this, over smoothness are often found in BMSR which tends to lose some important image details.

Evidences show that the combination of local sparsity and nonlocal self-similarity can be of potential power to further improve image quality for image restoration and recovery [14, 15]. The method of nonlocal means (NLM) filtering [16] has shown promising advantages in some typical tasks [17, 18] as well as in CT imaging applications [19, 20]. The NLM weights the pixels according to the degree of similarity among surrounding pixels which can be of well help to sharpen edges, suppress noise and preserve details for medical images. However, typical nonlocal methods, e.g., BMSR and NLM method results a non-convex functional which implies that only application of single nonlocal method may result unstable or not convergent algorithm. Nonetheless, the coupling of nonlocal regularization and local methods can be of a promising approach to improve reconstruction quality while with balance between stableness and performance. The work of Ertas [21, 22] proposed to used TV and NLM as the cascaded post-processing steps following the algebraic reconstruction technique (ART). However, this kind of combination does not have clearly defined optimization model, and also it is hard to obtain (nearly) theoretical results.

In this work, a combination of a typical representative local regularization method (TV) and a non-local method (NLM) is proposed, which is expected to be better than a purely local method or a purely non-local method. This combination is not only an improvement of TV, but also is regarded as a new method of combining local and non-local methods to design regularization. To the best of our knowledge, this is the first time that NLM combined with TV is applied in CT image reconstruction with clearly defined optimization model and efficient solvation algorithm developed by ADMM. In order to avoid to computing the inverse of a huge and dense matrix containing the system matrix, a first order linear expansion at the intermediate point is applied, together with Fast Fourier Transform (FFT) is utilized to tackle with the inverse of the composition of gradient operators.

The rest of this paper is organized as follows. Section 2 introduces the regularization of nonlocal means filtration and the optimization of total variation combining nonlocal filtration regularization, as well as an efficient algorithm developed by ADMM framework. Section 3 reports the experiments on simulation and real projection datasets, together with the comparisons with typical competing algorithms. Brief discussion and conclusion are presented in section 4 and section 5, respectively.

2 Method

2.1 Regularization by nonlocal filtration

We denote the nonlocal means filtering on an image x is as f_nlm (x), and for each pixel at i in x, the processing can be written as ${\hat{x}}_{i} = f_{nlm} (x_{i}),$ (1) where f_nlm (x_i) = ∑_{j∈T_j}w (i, j) x (j) and the weights w (i, j) | _{j∈T_j} is the similarity between the pixel i and j, satisfying 0 ≤ w (i, j) ≤1 and ∑_{j∈T_j}w (i, j) = 1, and T_k denotes a square search window of size of t_s × t_s centered at pixel k. The similarity w (i, j) is computed based on the pixels values of x (B_i) and x (B_j), where B_k denotes a square local patch with size of b_s × b_s centered at a pixel k, specifically, the corresponding weights w_ij defined as $w (i, j) = exp (- {∥ x (B_{j}) - x (B_{i}) ∥}_{2, σ}^{2} / h^{2}) / Z_{i},$ (2) where Z_i is the normalizing factor $Z_{i} = \sum_{j} exp (- {∥ k_{j} - k_{i} ∥}_{2, σ}^{2} / h^{2})$ , and ${∥ k_{j} - k_{i} ∥}_{2, σ}^{2} = \sum_{s} exp (- \frac{d^{2} (s)}{2 σ^{2}} {(k_{j} (s) - k_{i} (s))}^{2})$ is the Gaussian weighted Euclidean distance with zero mean and standard deviation of σ, and d (s) is pixel distance from the pixel s to central pixel x_i.

An equivalent matrix formation for (1) can be expressed as $\hat{x} = W_{σ} x, W_{σ} (i, j) = {\begin{matrix} w_{ij}, u_{j} \in k_{i} \\ 0, u_{j} \notin k_{i} \end{matrix} .$ (3)

In the following context, W_σ is simply denoted as W when not causing misunderstanding. With the above definitions and notations, a nonlocal means filtered image prior-based image regularization can be formulated as $R (x) = {∥ x - Wx ∥}_{2}^{2},$ (4) where the coefficient matrix W is calculated under given image x, search size of window, control factor h and the normalized factor c_i. The sizes of the local and search window have important impacts on the computation of the NLM-based regularization. Usually, large size of local and search window results relatively better image quality, however bringing heavy computation burden.

On one hand, the NLM-based regularization term R (x) has some properties similar to total variation on piece-wise constant image. For a sub-region r ⊆ x of an image x with all entities in r shares the same value, i.e., r_i = v for all i in this sub-region, it is easy to see that all weights are the same by (2) resulting Wx = x and R (r) =0. On the other hand, however, R (x) has some properties very different to that of total variation on non-piecewise constant images. For an image with uniform changing pixel values, the NLM filtering still holds Wx = x or Wx ≈ x which makes R (x) equals (or very close) to 0. However, this is not the case for total variation regularization.

2.2 Optimization by total variation combining nonlocal filtration regularization

According to the discrete-to-discrete (DD) model, the image to be reconstructed can be approximately represented by a vector as $x \in ℝ^{n}$ where n is the total number of pixels in the image. The measurements $b \in ℝ^{m}$ is obtained from the CT scanner where m is the number of total measurements. The scanning process is then approximately described by a linear transform of $A \in ℝ^{m \times n}$ as $Ax = b + r .$ (5)

Usually, the attenuation coefficients in $x \in ℝ^{n}$ are nonnegative real numbers and the elements in matrix A are generated using the intersection length of the ray and each entity in x by the Siddon’s model [23]. In the above equation, system matrix is known according to the hardware settings which is characterized by some critical parameters, i.e., distances from source to rotation center and detector, scanning views and positions and detector parameters. The projection dataset b is also known by observation. The core task of image reconstruction is to recover the intended image x with given A and b.

We proposed the total variation combining nonlocal means regularization as ${∥ Du ∥}_{1} + α {∥ u - Wu ∥}_{2}^{2},$ (6) subject to $Au - b = r, ∥ r ∥ \leq δ, u_{i} \geq 0 .$ (7)

With introducing new variables Du = z, u = x and penalize the constraints to the objective resulting $(z, u, x, λ, γ, σ) = \underset{z, u, x, λ, γ, σ}{arg min} {\begin{matrix} {∥ z ∥}_{1} + \frac{β}{2} {∥ z - Du - \frac{λ}{β} ∥}^{2} + \frac{μ}{2} {∥ Au - b - r - \frac{γ}{μ} ∥}^{2} \\ + α {∥ x - Wx ∥}_{2}^{2} + \frac{ρ}{2} {∥ u - x - \frac{σ}{ρ} ∥}^{2} \end{matrix}},$ (8) where λ, γ, σ are multipliers introduced by the augmented Lagrangian function. The updating of variables and multipliers can be expressed as $(z, u, x) = \underset{z, u, x}{arg min} {\begin{matrix} {∥ z ∥}_{1} + \frac{β}{2} {∥ z - Du - \frac{λ}{β} ∥}^{2} + \frac{μ}{2} {∥ Au - b - r - \frac{γ}{μ} ∥}^{2} \\ + α {∥ x - Wx ∥}_{2}^{2} + \frac{ρ}{2} {∥ u - x - \frac{σ}{ρ} ∥}^{2} \end{matrix}},$ (9) and ${\begin{matrix} λ \to λ - β (z - Du) \\ γ \to γ - μ (Au - b - r) \\ σ \to σ - ρ (u - x) \end{matrix} .$ (10)

Inspecting the above function (9), we find there are two blocks of variables, i.e., {z, x} and {u} which indicates the alternating direction method can be applied to solve the problem with guarantee of convergence.

2.3 Practical algorithm by ADMM

The problem of (9) can be solved by alternating (10) and solving the following three sub-problems as $\begin{matrix} \underset{z}{arg} min {∥ z ∥}_{1} + \frac{β}{2} {∥ z - Du - \frac{λ}{β} ∥}^{2}, \\ \underset{u}{arg} min \frac{β}{2} {∥ z - Du - \frac{λ}{β} ∥}^{2} + \frac{μ}{2} {∥ Au - b - r - \frac{γ}{μ} ∥}^{2} + \frac{ρ}{2} {∥ u - x - \frac{σ}{ρ} ∥}^{2}, \\ \underset{x}{arg} min α {∥ x - Wx ∥}_{2}^{2} + \frac{ρ}{2} {∥ u - x - \frac{σ}{ρ} ∥}^{2} . \end{matrix}$ (11)

The problem of z can be easily solved by the soft-shrinkage operator as $z \to max (| Du + \frac{λ}{β} | - \frac{1}{β}, 0) \cdot sign (Du + \frac{λ}{β}),$ (12) where sign is the signum function and the ^″ · ″ is the component wise operation.

For the sub-problem of u, the variable u are coupled with two operator of D and A. In this work, we introduce an efficient algorithm with can avoid the computation of the inverse of a huge matrix involving A^TA. We apply a first order approximation at the current point of u^k for the quadratic term involving Ax as $\frac{1}{2} {∥ Au - b - r - \frac{γ}{μ} ∥}^{2} \approx \frac{1}{2} {∥ A u^{k} - b - r - \frac{γ}{μ} ∥}^{2} + g_{k}^{T} (u - u^{k}) + \frac{1}{2 τ} {∥ u - u^{k} ∥}^{2},$ (13) where g_k = A^T (Au^k - b - r - γ/μ) is the gradient at u^k. The introduced parameter τ should be τ ≥ 1.0, and in this work it is set as τ = 1.0. Plug right hand side part of the above approximation formulation into (11) resulting $\underset{u}{arg} min \frac{β}{2} {∥ z - Du - \frac{λ}{β} ∥}^{2} + μ g_{k}^{T} (u - u^{k}) + \frac{μ}{2 τ} {∥ u - u^{k} ∥}^{2} + \frac{ρ}{2} {∥ u - x - \frac{σ}{ρ} ∥}^{2} .$ (14)

The first order optimal condition of the above function can be expressed as $β D^{T} (Du + \frac{λ}{β} - z) + μ A^{T} (A u^{k} - b - r - \frac{γ}{μ}) + \frac{μ}{τ} (u - u^{k}) + ρ (u - x - \frac{σ}{ρ}) = 0$ (15)

and also equivalently as $(β D^{T} D + \frac{μ}{τ} I + ρ I) u = rh s^{k},$ (16) where rhs^k = μu^k/τ - A^T (μAu^k - μb - μr - γ) + D^T (βz - λ) + ρx + σ. Under the periodic condition of the finite difference operator D, the above equation can be efficiently solved by FFT as $u = Re {F^{- 1} (F (rh s^{k}) . / F (β D^{T} D + μ I / τ + ρ I))},$ (17) where the division./ is an component wise operation. The operation of Re { c } returns the real part of a complex number c. Under the consideration of the constraints of u_i ≥ 0, all negative elements in u are forced to be 0.

For the problem of x, replace y = u - σ/ρ and this sub-problem can be equivalently expressed as $\underset{x}{arg} min \frac{α}{ρ} {∥ x - Wx ∥}_{2}^{2} + \frac{1}{2} {∥ x - y ∥}^{2} .$ (18)

The weight matrix W is calculated according to the input vector r in (18) and updated iteration by iteration using (2). And with the coefficient $\frac{α}{ρ} ⪡ 1.0$ , the dominant term in the above problem is $\frac{1}{2} {∥ x - y ∥}^{2}$ , which means the input y and the output x can be approximately deemed as very close to each other. Therefore, replace Wx with Wy, resulting $\underset{x}{arg} min \frac{α}{ρ} {∥ x - Wy ∥}_{2}^{2} + \frac{1}{2} {∥ x - y ∥}^{2} .$ (19) and the results of x can be written as $x = (2 α Wy + ρ y) / (2 α + ρ) .$ (20)

It must be pointed out that, the computation of nonlocal filtration of y (i.e, Wy) in (20) is actually of high complexity since its needs nonlocal searching on every pixels. In this work, a graphics processing unit (GPU) based parallel acceleration [24] is applied to improve the efficiency. The residual r between the estimated data and the observed data can be updated by the projection onto a L2-norm ball as $r^{k + 1} \to min (∥ Ax - b ∥, δ) \cdot \frac{(Ax - b)}{∥ A x^{k + 1} - b ∥},$ (21) where min(x₁, x₂) returns minimum value of the two scalars x₁, x₂. As for the sub-problems of (u, z, x) has already been solved with efficient methods, the proposed NLM-TV algorithm can be listed in Table 1 as follows. The stop condition is satisfied when the relative change ∥x^k+1 - x^k ∥/∥ x^k ∥ ≤ e_rel or the iteration reach the maximum number of loops.

Table 1

Pseudo-code of the proposed NLM-TV algorithm

Algorithm: nonlocal means filtration and total variation minimization

Input:A, b, β → β₀, μ → μ₀, α → 0.02, ρ → 1.0, z^k → 0,

u^k → 0, x^k → 0, r^k → 0, λ^k → 0, γ^k → 0, σ^k → 0, δ → δ₀, k → 0,

Repeat:

Update z: z^k+1 → max(|Du^k + λ^k/β| - 1/β, 0) · sign (Du^k + λ^k/β) ,

Update u: u^k+1 = Re{ F^-1 (F (rhs^k) ./F (βD^TD + μI/τ + ρI)) } with all negative elements forced to be 0.

Update weight matrix elements:

w_{ij} = exp (- {∥ k_{j} - k_{i} ∥}_{2, σ}^{2} / h^{2}) / Z_{i} .

Update x: x^k+1 = (2αWy^k + ρy^k)/(2α + ρ),

Update residual: r^k+1 → min(∥ Ax^k+1 - b ∥ , δ) · (Ax^k+1 - b)/∥ Ax^k+1 - b ∥ ,

Update multipliers:

λ^k+1 → λ^k - β (z^k+1 - Du^k+1) ,

γ^k+1 → γ^k - μ (Au^k+1 - b - r^k+1)

σ^k+1 → σ^k - ρ (u^k+1 - x^k+1),

k → k + 1,

When the stop condition satisfied, then stop and return u as the output x.

In the practical implementation of solving the sub-problem of x, there are some issues related to the properties of W in (20). It is obviously that the concrete form of W depends on the input of y which means that W will be changing during iterations resulting the problem of x is non-convex. However, there are two remedies to this concern. The first is to choose proper initialization of the starting point of iterations for u, and in this work u⁰ → 0 performs well. The second is to control the coefficient 0 < α/ρ ⪡ 1.0 for it is obvious that lim _α/ρ→0 (2αW + ρI)/(2α + ρ) = I. It must be pointed out that, in k-th iteration, the NLM filtering only applies on the input intermediate input (denoted as y^k). Therefore, the NLM filtering matrix may be denoted as W^(k). As the updating performs in an asymptotic fashion, it is reasonable to express W^k+1 as W^k+1 = W^k + Δ^k where Δ^k is a small perturbation of W^k+1 compared to W^k. With proper choices of starting point of iterations and controlling the coefficient 0 < α/ρ ⪡ 1.0 (α/ρ = 0.03 in this paper), it is reasonable to assume that the difference of y^k+1 over y^k is bounded, i.e., $∥ Δ y^{k} ∥ = ∥ y^{k + 1} - y^{k} ∥ \leq ɛ_{y}, ɛ \geq 0,$ (22) which implies that the perturbation Δ^k also bounded as ∥Δ^k ∥ ≤ ɛ_Δ 1 where ɛ_Δ ≥ 0, since W^k is determined by y^k.

Consider the difference of x^k+1 over x^k, i.e., Δx^k = x^k+1 - x^k whose modulus can be expressed as $∥ Δ x^{k} ∥ = ∥ α_{1} (W^{k} \cdot Δ y^{k} + Δ^{k} \cdot y^{k} + Δ^{k} \cdot Δ y^{k}) + α_{2} Δ y^{k} ∥,$ (23) where α₁ = 2α/(2α + ρ) , α₂ = ρ/(2α + ρ). Furthermore, the upper bound of ∥Δx^k∥ can be estimated as $∥ Δ x^{k} ∥ \leq α_{1} ∥ W^{k} \cdot Δ y^{k} ∥ + α_{1} ∥ Δ^{k} \cdot y^{k} ∥ + α_{1} ∥ Δ^{k} \cdot Δ y^{k} ∥ + α_{2} ∥ Δ y^{k} ∥,$ (24) and furthermore resulting $∥ Δ x^{k} ∥ \leq α_{1} ɛ_{y} ∥ W^{k} ∥ + α_{1} ɛ_{Δ} ∥ y^{k} ∥ + α_{1} ɛ_{y} ɛ_{Δ} + α_{2} ɛ_{y},$ (25) where elements in W^k satisfying 0 ≤ w (i, j) ≤1 and ∑_{j∈T_j}w (i, j) = 1 obviously resulting W^k a bounded matrix. Therefore, with (25), the difference of x^k+1 over x^k is also bounded. This implies that the iteration with the NLM filtration is relatively stable, which is actually also verified by the experiments in the following section.

3 Experiments and results

3.1 Tests on simulation dataset

In this section, a series of simulation experiments are carried out on a chest phantom as shown in Fig. 2 (left). The grey value of the chosen phantom is scaled in the range of [0.0, 1.0]. Filtered back-projection (FBP) reconstruction and four typical iterative algorithms, i.e., simultaneous algebra reconstruction technique (SART), SART with total variation regularization (SART-TV), TV minimization by ADMM (TV-L1) and the recently proposed nonlocal block-matching sparsity regularization (BMSR), are included as the competing methods to the proposed method in this work.

Fig. 2

Reference image (left) and the corresponding 19-view FBP reconstruction (right) with display grey scale in [0.1 0.5].

All the experiments are conducted on a laptop configured with a dual-core CPU (@ 1.60 GHz, only single core used) with memory of 8.00 GB. An Nvidia GPU (TITAN V with compute capability of 7.0, total memory of 12.00 GB and 5120 CUDA cores) is utilized to accelerate the NLM filtration (20) for improving the computation efficiency. All the initial guesses of the iterative methods are set to zero.

3.1.1 Parameter settings

In the settings of the simulation tests on the chest phantom, the scanning geometry are circular trajectory. The distances from the X-ray source to the rotation center and the detection center are 503.1 mm and 1052.8 mm. The scanning view number is 19 evenly distributed within 2π. The number of detector bins is 512 with size of 0.628 mm for each of them. The size of reconstruction image is 256 by 256, with each pixel of 0.60 mm by 0.60 mm.

As for the settings of the proposed algorithm, the typical parameters for the proposed NLM-TV are α, β, μ and ρ. It should be noted that the optimal settings for parameters is a non-trivial work. However, the empirical settings for these are chosen as β = 32, μ = 1024 and α/ρ = 0.03 in this work. The choice of δ₀ is actually determined by the noise level of the observation projection data. Therefore, in the simulation test, δ₀ is set to 0 since we only focus on sparse-view reconstruction under a noiseless simulation.

For the parameters (b_s, t_s, σ, h) of NLM, for images with the dynamic range of [0.0 1.0], an empirical and proper setting for σ, h is chosen and fixed as σ = 1.0 and h = 0.02. As for local patch size b_s and the nonlocal search window size t_s, a serial tests are carried out on the simulation tests with different pares of b₃ and t₃. A series of plots of RMSEs on different b₃ and t₃ are drawn in Fig. 1 in which a basic conclusions can be that, with larger b₃ and t₃, the reconstructions provide better RMSEs as the functions of iteration loops (left of Fig. 1), but relatively slower RMSEs as the functions of time consumption (right of Fig. 1). For the comprehensive capabilities of the TV-NLM method, the pair of b₃ = 7 and t₃ = 13 is chosen and fixed in the following experiments.

Fig. 1

RMSEs of different patch size and search window size as functions of iteration loops (left) and time consumption (right) of the proposed TV-NLM method.

The parameters for the chosen competing methods are tuned to be optimized for a relatively fair comparison. For simulation studies, the relative change tolerance e_rel is set to 10^-5. And the maximum iteration number is set to 40,000 to give a relatively full investigations of the iterative behaviors for all the iterative algorithms. The running time for each iterative algorithm are another significant index. According to the sampling condition for the simulation dataset, the average computation for SART, SART-TV, TV-L1, BMSR and the propose TV-NLM are about 0.03 seconds, 0.37 seconds, 0.025 seconds, 0.075 seconds and 0.06 seconds for a single iteration loop, respectively.

3.1.2 Results of simulation studies

The analytical by FBP (with Shepp-Logan filter) reconstruction are shown in right of Fig. 2, which is severely degraded by the streak artifacts due to the data deficiency. The background truth and the reconstructions by iterative algorithms are shown in Fig. 3 displayed in gray window of [0.3 0.48]. The result of SART is still of very low quality showing lots of noise and artifacts. The results under regularizations (c–f of Fig. 3) of TV, block matching and nonlocal means are of better quality. Among the iterative results, regularizations with only TV suffers from stair-stepping artifacts. The result of the BMSR shows better performance on denoising while has drawbacks of over smoothing, which sacrifice much detail structures. The reconstruction of TV-NLM shown in Fig. 3 (f) has comprehensively advantages on noise suppression and detail preservation over SART-TV, TV-L1 and BMSR. Since there are much details in the lung tissue, therefore, a region-of-interest of the lung area which is marked by a rectangle in Fig. 2 for further visual inspection. A low and narrower intensity gray window is chosen to display the results of iterative algorithms, shown in Fig. 4. It is obvious that the proposed method shows clearer details structures than those by BMSR. Moreover, the error maps x - x₀ (where x is the reconstruction and is the background truth) for each iterative algorithm are displayed in Fig. 5.

Fig. 3

Reconstructions of iterative algorithms under 19-view projections. Images from (a)-(f) are background truth, and results of SART (b), SART-TV (c), TV-L1 (d), BMSR (e) and the proposed TV-NLM (f), respectively. The displayed gray windows are all [0.3, 0.48].

Fig. 4

Reconstructions of iterative algorithms under 19-view projections. The same as Fig. 2 but displayed in gray window of [0.01, 0.2].

Fig. 5

Error map of each iterative algorithm under 19-view projections. The same as Fig. 2 but displayed in gray window of [–0.3, 0.3].

In order to numerically compare the results, root mean squared error (RMSE), structure similarity of image (SSIM) and the peak signal-to-noise ratio (PSNR) are calculated. The smaller the RMSE value or the larger the PSNR value, the smaller the distortion of the reconstructed image. Moreover, the SSIM indicator measures image similarity based on three features: brightness, contrast, and structure. Generally, SSIM is in the range [0, 1]. The larger the value, the smaller the image distortion. The optimal values of RMSE, SSIM and PSNR during the iterations of each algorithm are listed in Table 2. It should be noted that the proposed new method shows obvious advantages over the competing methods in terms of RMSE, SSIM and PSNR. For the convergence comparison of the iterative behaviors, RMSEs of intermediate results of typical iterative number are calculated and listed in Table 3. Moreover, the RMSEs of each algorithm under iteration number are plotted in Fig. 6. The numerical results (by Table 2, Table 3 and Fig. 6) verify that TV-NLM algorithm outperforms the competing method under the reconstruction from 19-view projections.

Table 2

Optimal numerical criteria for iterative reconstructions

Index	SART	SART-TV	TV-L1	BMSR	TV-NLM
RMSE	0.0338	0.0182	0.0133	0.0125	0.0093
SSIM	0.6386	0.08765	0.9170	0.9264	0.9528
PSNR	26.3103	31.6914	34.1825	34.9458	36.9216

Table 3

Iteration behaviors for iterative algorithms

RMSE	100	1000	5000	10000	20000	30000	40000	optimal
SART	0.0346	0.0340	0.0338	0.0338	0.0339	0.0342	0.0346	0.0338
SART-TV	0.0203	0.0182	0.0182	0.0182	0.0182	0.0182	0.0182	0.0182
TV-L1	0.0219	0.0134	0.0133	0.0133	0.0133	0.0133	0.0133	0.0133
BMSR	0.0323	0.0229	0.0200	0.0178	0.0126	0.0125	0.0125	0.0125
TV-NLM	0.0175	0.0095	0.0096	0.0096	0.0096	0.0096	0.0097	0.0093

Fig. 6

Iterative comparison of iterative algorithms for reconstructions from 19-view projection.

3.2 Verifications on realistic dataset

In order to verify the behaviors of the proposed on real dataset, our laboratory uses an industrial CT system to conduct experiments on real dataset. The industrial CT system in our laboratory is a CT system for circular trajectory scanning, which includes detectors, ray sources, stages, mechanical systems, and computer acquisition and projection software. The ray source is 150 kV MICROFOCUS X-RAY SOURCE L12161-07, and the detector is THALES PIXIUM RF 4343. First, we put the dead rabbit head into a plastic container for experimentation. It scans 360° along a circular trajectory with a tube voltage of 80 kVp and a tube current of 125μA. The distance from the source to the rotation center is 502.8 mm, and the distance from the source to the center of the detector is 1434.7 mm. Only the central slice of each cone beam projection is taken into account for 2D reconstruction, thus making the scanning to be the fan beam geometry. There are 768 bins for each view, and the equivalent pixel size of each bin is 0.508 mm. A total of 360 views within 360 degrees are acquired and SART is applied on 360 views to obtain a reference reconstruction. For the sparse view reconstruction, 36 and 20 projections are evenly chosen with angular step of 10 and 18 degrees within 360 degrees, respectively.

3.2.1 Results of 36 views of real projections

The reconstruction size of images are all 384 by 384, with each pixel of 0.3561 mm by 0.3561 mm. Under real projection dataset of 36 views, the average computation for SART, SART-TV, TV-L1, BMSR and the propose TV-NLM are about 0.08 seconds, 0.84 seconds, 0.04 seconds, 018 seconds and 0.07 seconds per iteration loop, respectively. It must be pointed out that the real CT projections contain noise or inconsistencies in the observation which cannot be avoid. Therefore, the iteration is terminated earlier than simulation dataset. The relative tolerance e_rel is set to 10^-4 and the maximum iteration number is set to 2,000 for all iterative algorithms.

The real data results of 36 views of projections are displayed in Fig. 7, where the error maps of all the iterative reconstructions are shown in Fig. 8. The results in Fig. 7 indicates that TV-based methods (SART-TV and TV-L1, i.e., Fig. 7-c and Fig. 7-d) provide much better outcomes than algebraic algorithm (Fig. 7-b). Although the streak artifacts are efficiently suppressed by the TV methods, yet TV suffers from staircase artifacts. The results of BMSR methods (Fig. 7-e) improves the staircase artifacts compared with TV-base ones. However, the over smoothing are presented in Fig. 7-e, which sacrifices lots detail structures compared with the reference provided by SART under 360-view projections. The results provided by the proposed TV-NLM method in Fig. 7-f balanced with detail structure preservation the staircase artifacts suppression.

Fig. 7

Results of reference and 36-view reconstructions: (a) SART reference (with 360 views of projections), (b–f) results of SART (b), SART-TV (c), TV-L1 (d), BMSR (e) and the proposed TV-NLM (f) with 36 views of projections, respectively. The displayed windows are all [0.002, 0.06] mm^–1.

Fig. 8

Error maps of iterative algorithms, the same layout arrangement as Fig. 7 but in a display window of [–0.025 0.025].

The error maps presented in Fig. 8 also give an intuitive sense of the relative error of each method compared to the reference. In order to make a further inspection of the results, a region-of-interest (ROI) is chosen and its position is marked by a red solid rectangle in Fig. 7 and their enlarged version of each method are displayed in Fig. 9. From the visual inspection of these enlarged ROIs, a basic property of the proposed method can be drawn that the TV-NLM gives better balance between detail preservation the noise suppression. For quantitive comparisons, three typical indices including RMSE, PSNR and SSIM are calculated for full images and the selected ROIs and listed in Table 4.

Fig. 9

Enlarged display of real data results (reference and 36-view iterative reconstructions) in the selected area (marked by rectangle in Fig. 7). The displayed windows are all [0.002, 0.06] mm^–1.

Table 4

Numerical criteria for iterative reconstructions under real projection dataset of 36 views

Area	Index	SART	SART-TV	TV-L1	BMSR	TV-NLM
Full	RMSE	0.0049	0.0032	0.0022	0.0022	0.0019
	PSNR	32.2523	35.9349	39.2866	39.1063	40.5783
	SSIM	0.8461	0.9537	0.9618	0.9643	0.9768
ROI	RMSE	0.0084	0.0060	0.0038	0.0035	0.0032
	PSNR	27.5589	30.3976	34.3463	35.0435	36.0504
	SSIM	0.7452	0.9059	0.9445	0.9514	0.9624

3.2.2 Results of 20 views of real projections

In order to test the proposed method with more ill-posed imaging situation, the comparisons under 20 views of real projection dataset are also carried. The reconstructions are displayed in Fig. 10 with the error maps shown in Fig. 11, where it is obvious that the results from 20 views are of lower quality than those from 36 views. Nevertheless, the result of TV-NLM method shows advantage of staircase artifacts suppression over TV-based methods in the enlarged ROIs shown in Fig. 12. Also, the compared with the BMSR in Fig. 12-e, the proposed method provides closer results to the reference. With further inspection into Fig. 12, it is quite obvious that BMSR method has strong capability in noise suppression, however, over smoothing erases some important details. The quantities for evaluating the image quality on full images and the selected ROIs are listed in Table 5. Both the visual inspections and the quantity evaluation indicate the advantages of recovery capability over the competing methods.

Fig. 10

Results of reference and 20-view reconstructions: (a) SART reference (with 360 views of projections), (b–f) results of SART (b), SART-TV (c), TV-L1 (d), BMSR (e) and the proposed TV-NLM (f) with 20 views of projections, respectively. The displayed windows are all [0.002, 0.06] mm^–1.

Fig. 11

Error map of iterative algorithms, the same as Fig. 10 but in a display window of [–0.025 0.025].

Fig. 12

Enlarged display of real data results (reference and 20-view iterative reconstructions) in the selected area (marked by rectangle in Fig. 7). The displayed windows are all [0.002, 0.06] mm^–1.

Table 5

Numerical criteria for iterative reconstructions under real projection dataset of 20 views

Area	Index	SART	SART-TV	TV-L1	BMSR	TV-NLM
Full	RMSE	0.0064	0.0042	0.0035	0.0032	0.0029
	PSNR	29.9269	33.5494	35.2265	35.8064	36.8411
	SSIM	0.7984	0.9208	0.9324	0.9486	0.9593
ROI	RMSE	0.0102	0.0066	0.0051	0.0043	0.0041
	PSNR	25.8336	29.6212	31.8072	33.3383	33.6694
	SSIM	0.7040	0.8740	0.9157	0.9279	0.9465

4 Discussion

In this work, a regularization based on total variation combing non-local means filtration has been proposed. It is obvious that the total variation is a kind of local operator-based method in which it only concerns the adjacent two pixels while the non-local means based method enlarged image pixels to a wider range. In practice, the total variation based methods are easily to generate stair-case or blocky artifacts. On the contrary, the block matching based method proposed recent years is a typical non-local regularization, however, it is of relatively high computation complexity and it also tends to over smoothing. In this work, the efforts of attempt of combining the two kinds of regularization have been made. The design of the proposed method is inspired by the combination of local operation and non-local operation aiming at providing improved results over the usage of single approach.

In this study, the ADMM framework is tailored to devise a practical algorithm for the proposed regularized optimization. With careful inspection of the regularization problem, the three variables can be actually divided into two blocks, i.e., (z, x) and u, which is strictly falls into the scope of ADMM with guarantee of convergence. As for the solvation of sub-problems with respect to each variable, soft shrinkage operation on z, first order linear expansion at current point of u and approximated transformation on x are applied with resulting an efficient algorithm. A series of experiments containing tuning of key NLM parameters, experiments of simulation dataset and verifications on real projection dataset are carried out. All the results indicate that the proposed method is of promising utility for sparse view CT image reconstructions compared to the competing methods, especially on noise suppression and detail preservation.

Moreover, an open source implementation of the proposed method will be provided by the authors for the interested readers and can be accessed online 2 in the future, with a demonstration on the simulation dataset studied in this paper. While the presented work is studied in two dimensional setting, however, it is quite conceptually, mathematically and practically straightforward that the proposed method can be directly extended into three dimension (3D) reconstruction. Furthermore, the 3D versions for the difference operator D and the NLM filtration, as well as the interface for generating 3D cone-beam system matrix have also been implemented and included in the open source implementation. Interested readers are encouraged to conduct more investigations.

5 Conclusion

By combing the total variation and the non-local means filtration, in this study, a novel image reconstruction method has been proposed. The alternating direction method of multipliers is applied to devise a practical algorithm, with efficient solvers for each sub-problem has been developed. A series of experiments, including simulation and real projections dataset, have been conducted, and the comparisons with the selected competing methods verifies the capabilities of the proposed method. The investigations also verify the advantages of the novel method. The authors are going to conduct further investigations of the proposed method and extend to cone-beam scanning setting for more realistic applications.

Footnotes

Acknowledgments

This study is supported by the National Natural Science Foundation of China (NSFC), No. 62101596 and the China Postdoctoral Science Foundation (CPSF), No. 2019M663996. This work is also supported by the National Key Research and Development (NKRD) Project of China, No. 2020YFC1522002. The authors claim that the points of views in this work are solely responsible to the authors themselves and not necessarily the official views of the NSFC, CPSF and NKRD. The authors would like to thank the anonymous reviewers for their helpful and constructive suggestions.

Competing interests

The authors declare that they have no competing interests.

Matrix norm of A induced by vector norm of x is defined as ∥A ∥ = sup {∥ Ax ∥ : ∥ x ∥ = 1} or equivalently as ∥A ∥ = sup {∥ Ax ∥/∥ x ∥ : x ≠ 0}.

References

Geyer

L.L.

, Schoepf

U.J.

, Meinel

F.G.

, et al., State of the art: iterative CT reconstruction techniques, Radiology 276(2) (2015), 339–357.

Jorgensen

J.S.

, Sidky

E.Y.

and Pan

, Quantifying admissible undersampling for sparsity-exploiting iterative image reconstruction in X-ray CT, IEEE Transactions on Medical Imaging 32(2) (2013), 460–473.

Engl

H.W.

, Hanke

and Neubauer

, Regularization of inverse problems, Springer Science & Business Media, 1996.

Donoho

D.L.

, Compressed Sensing, IEEE Transactions on Information Theory 52(4) (2006), 1289–1306.

Sidky

E.Y.

and Pan

X.C.

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

Qiao

and Lu

, A TV-minimization image-reconstruction algorithm without system matrix, Journal of X-Ray Science and Technology 29(5) (2021), 851–865.

Qiao

, A simple and fast ASD-POCS algorithm for image reconstruction, Journal of X-Ray Science and Technology 29(3) (2021), 491–506.

Wahlberg

, Boyd

, Annergren

, et al., An ADMM algorithm for a class of total variation regularized estimation problems. arXiv: 12031828, 2012,

Zhang

H.M.

, Wang

L.Y.

, Yan

, et al., Image reconstruction based on total-variation minimization and alternating direction method in linear scan computed tomography, Chinese B. Physics 22(7) (2013), 582–589.

10.

Kim

, Chen

, Wang

, et al., Non-local total-variation (NLTV) minimization combined with reweighted L1-norm for compressed sensing CT reconstruction, Physics in Medicine and Biology 61(18) (2016), 6878–6891.

11.

Andersen

and Kak

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

12.

Cai

, Li

, Zheng

, et al., Block matching sparsity regularization-based image reconstruction for incomplete projection data in computed tomography, Physics in Medicine & Biology 63(3) (2018), 035045.

13.

Cai

, Li

, Zheng

, et al., Block-matching sparsity regularization-based image reconstruction for low-dose computed tomography, Medical Physics 45(6) (2018), 2439–2452.

14.

Ertas

, Akan

, Yildirim

, et al., Image denoising by using non-local means and total variation, 307 Proceedings of the 2014 22nd Signal Processing and Communications Applications Conference (SIU), F 23-25 April 2014, 2014.

15.

Zhang

, Liu

, Xiong

, et al., Improved total variation based image compressive sensing recovery by nonlocal regularization, Proceedings of the 2013 IEEE International Symposium on Circuits and Systems (ISCAS), F 19-23 May 2013, 2013.

16.

Buades

, Coll

and Morel

J.M.

, A review of image denoising algorithms, with a new one, Siam Journal on Multiscale Modeling & Simulation 4(2) (2005), 490–530.

17.

Lou

, Zhang

, Osher

, et al., Image recovery via nonlocal operators, Journal of Scientific Computing 42(2) (2010), 185–197.

18.

Lebrun

, Buades

and Morel

J.M.

, A nonlocal Bayesian image denoising algorithm, SIAM Journal on Imaging Sciences 6(3) (2013), 1665–1688.

19.

Zeng

, Huang

, Zhang

, et al., Spectral CT image restoration via an average image-induced nonlocal means filter, IEEE Transactions on Biomedical Engineering 63(5) (2016), 1044–1057.

20.

Zhang

, Zeng

, Zhang

, et al., Applications of nonlocal means algorithm in low-dose X-ray CT image processing and reconstruction: A review, Medical Physics 44(3) (2017), 1168–1185.

21.

Ertas

, Yildirim

, Kamasak

, et al., An iterative tomosynthesis reconstruction using total variation combined with non-local means filtering, Biomedical Engineering Online 13 (2014), 65.

22.

Ertas

, Yildirim

, Kamasak

, et al., Iterative image reconstruction using non-local means with total variation from insufficient projection data, Journal of X-Ray Science and Technology 24(1) (2016), 1–8.

23.

Siddon

R.L.

, Fast calculation of the exact radiological path for a three-dimensional CT array, Medical Physics 12(2) (1985), 252–255.

24.

Upadhya

A.H.

, Talawar

and Rajan

, GPU implementation of non-local maximum likelihood estimation method for denoising magnetic resonance images, Journal of Real-Time Image Processing 13(1) (2017), 181–192.