Combine TV-L1 model with guided image filtering for wide and faint ring artifacts correction of in-line x-ray phase contrast computed tomography

Abstract

In practice, mis-calibrated detector pixels give rise to wide and faint ring artifacts in the reconstruction image of the In-line phase-contrast computed tomography (IL-PC-CT). Ring artifacts correction is essential in IL-PC-CT. In this study, a novel method of wide and faint ring artifacts correction was presented based on combining TV-L1 model with guided image filtering (GIF) in the reconstruction image domain. The new correction method includes two main steps namely, the GIF step and the TV-L1 step. To validate the performance of this method, simulation data and real experimental synchrotron data are provided. The results demonstrate that TV-L1 model with GIF step can effectively correct the wide and faint ring artifacts for IL-PC-CT.

Keywords

Phase-contrast computed tomography guided image filtering TV-L1 model ring artifacts

1 Introduction

Phase-contrast imaging (PCI) is a new imaging method that can reconstruct proper image for soft tissue [1]. Additionally, PCI is considered a powerful preclinical imaging modality to observe fine structures with a resolution higher than any available clinical radiography [2]. Computed tomography is effective and efficient for noninvasive diagnosis and non-destructive technique. In recent years, by using PCI together with CT, the Phase contrast CT (PC-CT) technique is also developed, which possesses excellent advantages compared to the conventional absorption-based CT when imaging soft tissues [3 –5]. PC-CT methods include interferometric CT [6], diffraction-enhanced CT [7, 8], grating-based phase contrast imaging methods [9, 10], and in-line PC-CT (IL-PC-CT) methods [11, 12].

IL-PC-CT images, reconstructed from the 2D X-ray projection data acquired with a flat-panel detector, are often corrupted by ring artifacts. As visualization and quantification of the anatomic and pathological features are influenced by these ring artifacts, removal or reduction of these artifacts is essential. Ring artifact removal or reduction without impairing the image quality is still a challenging problem for the researchers. To overcome the challenge, several methods have been reported so far to correct the ring artifacts from CT images.

The flat field correction is a common method to reduce ring artifacts [13]. However, ring artifacts are not removed completely by this method when different camera elements have intensity dependent, nonlinear response functions or the incident beam has time dependent non-uniformities [1]. Another method for ring artifacts suppression is based on moving the sample or the detector system during the acquisition in defined horizontal and vertical steps [14, 15]. However, this requires very precise translation motor components. The third method for correct ring artifacts is the image-processing based on sinogram preprocessing and reconstructed images post-processing. Some reported methods dealt with the raw sinogram aim at detecting and correcting the spurious lines in the sinogram before applying the reconstruction process [1, 16]. The post-processing approach may be tempting to detect and remove the ring artifact in the image domain [17] rather than the sinogram domain [18 –20]. Differing from classical sinogram preprocessing and image postprocessing techniques, a novel approach is presented to correct ring artifacts in compressed sensing tomographic reconstruction [21]. It is shown that the ring artifacts correction in compressed sensing tomographic reconstruction can be easily integrated in the formalism, enabling simultaneous slice reconstruction and ring artifacts correction. The optimization model of compressed sensing tomographic reconstruction is solved by the fast iterative shrinkage-thresholding algorithms (FISTA).

The ring artifacts can be of different types and of different intensities. On the one hand, completely damaged detector pixels cause strong ring artifacts in the tomographic image [16]. On the other hand, mis-calibrated detector pixels, e.g. due to beam instabilities not completely taken into account by a normalization correction, give rise to wide and faint rings in the tomographic image [22]. The sinogram pre-processing technique is the wavelet-FFT filter [16]. Experiment result in [16] showed that most wide and faint ring artifacts were reduced from the reconstructed slice. However, the experiment results corrected by wavelet-FFT filter showed that there were minimal distortions at the reconstructed image center. In our study, we mainly research the correction method for wide and faint ring artifacts in the reconstruction image domain.

Recently, a guided image filtering (GIF) method was proposed that yielded both quality and efficiency in many applications [23, 24]. The guided filter has good edge-preserving smoothing properties, if a pixel is in the middle of a “high variance” area, then its value is unchanged, whereas if it is in the middle of a “flat patch” area, its value becomes the average of the pixels nearby [23]. In our study, the ring artifacts is characterized with wide and faint, so the wide and faint ring artifacts are in the “flat patch” area, the component of ring artifacts becomes the average of the pixels nearby. Thus, the GIF step can suppress the ring artifacts to some extent, however, the GIF step cannot removed the ring artifacts completely (see Section 3.1).

Sparse regularization techniques promote sparse solutions when solving ill-posed inverse problems. For example, L1 regularization tended to produce a sparse solution with clean background [25, 26]. However, in case of high fraction of noise, L1 regularization fails to recover the image efficiently. Therefore, we propose to add another TV regularization to correct ring artifacts. The TV term (isotropic or anisotropic) is essentially the L1-norm of the gradient of the image [27, 28], which can promote sparse gradients. The main idea of anisotropic TV is to separate the first-order finite difference operator in the horizontal and vertical directions to avoid the interference between them, hence the anisotropic total variation regularization can preserve image sharp with corners [29]. In industrial and clinical tomographic imaging applications, natural images usually contain a lot of anisotropic features and have a sparse representation in a transform domain. Therefore, we study the anisotropic TV regularization term added to L1-TV model.

Since we know the natural image is sparse in a transform domain and we expect the ring artifacts to not be sparse in the same transform domain (since it is not a real image). TV-L1 model not only can be viewed as promoting a doubly sparse structure respectively, but also combines the merits from both TV and L1 regularization in the sense that it successfully reconstructs piecewise constant image that preserve both geometry and details with little background noise [30, 31]. Thus, the TV-L1 model sparsity constraint forces the ring artifacts variables to have only a few not null components. Similarly, in this study, the ring artifacts can be of wide and faint, which can be corrected by TV-L1 model. Inspired by the work [23 , 31], the new ring artifact correction method is presented for IL-PC-CT based on combining TV-L1 model with guided image filtering (GIF). This new algorithm is referred to as the GIF-TV-L1 algorithm. GIF-TV-L1 algorithm can reduce visible wide and faint ring artifacts while preserve the image section details compared with wavelet-FFT filter. The proposed ring artifacts correction method works quite well in practice (see Section 3.2).

2 Method

2.1 Synchrotron IL-PC-CT data acquisition

This study was approved by the ethics committee of Tianjin Medical University, Tianjin, China, and written informed consent was obtained from all patients. The real liver sample experiments were performed in BL13W1 of the Shanghai synchrotron radiation facility (SSRF). Figure 1 shows the schematic of the IL-PC X-ray imaging setup. The monochromator crystal was a combination of a Si (111) orientation crystal and a Si (311) orientation crystal. A high-precision sample platform was used to position the sample and rotate the sample axis perpendicular to the beam. The x-ray beam energy in the experiments was set at 20 keV, and the detector employed an x-ray CCD camera system with 9×9μm² per pixel. The projections were recorded with a sample-to-detector distance of 1 m and exposure time of 10 ms. Projection data was collected from 1024 projections evenly distributed over pi-view in the BL13W1 of the SSRF. In addition, 20 flat field images and 10 dark field images were collected. The liver sample was placed on the shelf, and its image f^FBP of size 1024×1024 pixels was reconstructed using a FBP algorithm [13].

Fig.1

Schematic diagram of the in-line phase contrast imaging setup.

Figure 2(a) shows sinogram image of the projection data from the liver section. The data were collected by taking 1024 projections with total 180°angle of view. In Fig. 2(b), the magnified region of interest (ROI) of sinogram reveals vertical lines that are responsible for ring artifacts. In this study, our aim is to correct ring artifacts.

Fig.2

(a) Sinogram of the projection data. (b) magnified ROI in the Sinogram.

2.2 FBP-GIF algorithm

The guided image filtering [23] is a local affine model that connects the filtering output f^out to the guidance image f^guide: $f_{j}^{out} = a_{k} f_{j}^{guide} + b_{k}, \forall j \in w_{k}$ (1) where (a_k, b_k) are linear coefficients that are assumed to be constant in a square window w_k, r denotes a radius of a disc containing w_k. The values of a_k and b_k are then obtained by minimizing a cost function: $min_{f_{j}^{out}} \sum_{j \in w_{k}} ((a_{k} f_{j}^{guide} + b_{k} - f_{j}^{in})^{2} + ɛ a_{k}^{2})$ (2)

Where ɛ is a regularization parameter penalizing large a_k, fⁱⁿ is a filtering input image.

Equation (2) is the linear ridge regression model and its solution can be expressed by a linear regression [32]: $a_{k} = \frac{\frac{1}{| w |} \sum_{j \in w_{k}} f_{j}^{guide} f_{j}^{in} - η_{k} {\tilde{f}}_{k}}{σ_{k}^{2} + ɛ}$ (3) $b_{k} = {\tilde{f}}_{k} - a_{k} η_{k}$ (4)

Where η_k and $σ_{k}^{2}$ are, respectively, the mean and variance of f^guide in w_k, |w| is the number of pixels in w_k, and ${\tilde{f}}_{k} = \frac{1}{| w |} \sum_{j \in w_{k}} f_{j}^{in}$ is the mean of fⁱⁿ in w_k.

Because several windows w_k overlap a given pixel, so the value of $f_{j}^{out}$ in Equation (1) is not identical when it is computed in different windows. A simple strategy is to average all the possible values of $f_{j}^{out}$ . So after computing (a_k, b_k) for all the windows w_k in the image, the filter output is computed as ${f_{j}^{out}}_{=} \frac{1}{| w |} \sum_{k : j \in w_{k}} (a_{k} f_{j}^{guide} + b_{k}) = {\bar{a}}_{j} f_{j}^{guide} + {\bar{b}}_{j}$ (5) where ${\bar{a}}_{j} = \frac{1}{| w |} \sum_{k \in w_{j}} a_{k}$ and ${\bar{b}}_{j} = \frac{1}{| w |} \sum_{k \in w_{j}} b_{k}$ , $({\bar{a}}_{j}, {\bar{b}}_{j})$ are the average coefficients calculated from all the overlapping windows j.

Equations (3), (4), and (5) taken together define the guided image filtering, for which the main computational burden is a boxfilter of window radius r, denoted boxfilter. The guided image filtering is as follows: $f^{out} = boxfilter (a) \cdot * f^{guide} + boxfilter (b)$ (6) where ·* denotes element-by-element multiplication.

In our study, f^guide = fⁱⁿ, fⁱⁿ was reconstructed using a filtered back-projection (FBP) algorithm, this guided image filtering(GIF)algorithm is referred to as the FBP-GIF algorithm.

2.3 FBP-TV-L1 algorithm

Proper utilization of the prior information enables the image quality to be improved, sparse prior information has attracted considerable interest. The choice among the various prior terms depends on the specific object and on the desire to preserve or emphasize particular features [33, 34]. In order to combine the merits from both L1and TV regularization, the ring artifacts correction technique based on a combination of L1 and TV regularization is presented in IL-PC-CT. It is known that CT images usually have sparse representations under certain wavelet bases [35]. The Daubechies wavelet basis is used in our experiments. Let Ψ = [ψ₁, ψ₂,. . . . . . ψ_N] be the Daubechies wavelet basis. The ring artifacts correction technique becomes the TV-L1model: $min_{f} \frac{1}{2} {∥ f - g ∥}_{2}^{2} + λ_{1} {∥ f ∥}_{TV} + λ_{2} {∥ Ψ^{T} f ∥}_{1}$ (7)

For the sake of convenience, let ∥f ∥ _TV = ∥ Df ∥ ₁, then we can rewrite the sub-problem Equation (7) as $min_{f} \frac{1}{2} {∥ f - g ∥}_{2}^{2} + λ_{1} {∥ Df ∥}_{1} + λ_{2} {∥ Ψ^{T} f ∥}_{1}$ (8) where ∥Df ∥ ₁ denotes the anisotropic TV term; g is an input image corrupted by ring artifacts; λ₁ and λ₂ are the regularization parameters; ∥Ψ^Tf ∥ ₁ is the L₁-norm of the representation of f under Ψ. From the above equation it is obvious that TV-L1 model promotes simultaneous sparsity of the transform coefficients in the wavelet domain and the sparsity of the gradients in the TV domain. Hence, the TV-L1 model can be viewed as promoting a doubly sparse structure respectively.

ADMM is an algorithm that is intended to blend the decomposability of dual ascent with the superior convergence properties of the method of multipliers [36]. ADMM decomposes a large global problem into a series of smaller local sub-problems, and coordinates the local solutions to compute the globally optimal solution. Inspecting (8), we adopt a novel splitting strategy [31] by introducing three auxiliary variables y = f, Ψ^Tf = u, Df = z, the optimization is rewritten as $\begin{matrix} min_{f, y, z, u} \frac{1}{2} {∥ y - g ∥}_{2}^{2} + λ_{1} {∥ z ∥}_{1} + λ_{2} {∥ u ∥}_{1} \\ subject to y = f, z = Df, u = Ψ^{T} f \end{matrix}$ (9)

Its corresponding augmented Lagrangian function is $\begin{matrix} L (f, y, z, u, μ, ρ, γ) = \frac{1}{2} {∥ y - g ∥}_{2}^{2} + λ_{1} {∥ z ∥}_{1} + λ_{2} {∥ u ∥}_{1} \\ - μ^{T} (y - f) - ρ^{T} (z - Df) - γ^{T} (u - Ψ^{T} f) \\ + \frac{δ}{2} {∥ y - f ∥}_{2}^{2} + \frac{θ}{2} {∥ z - Df ∥}_{2}^{2} + \frac{η}{2} {∥ u - Ψ^{T} f ∥}_{2}^{2} \end{matrix}$ (10)

Where λ₁ and λ₂ are the regularization parameters; the vectors μ, ρ and γ are the Lagrange multipliers, and δ, θ and η are the penalty parameters.

We use the alternating direction method to iteratively solve the following sub-problems: $\begin{matrix} f^{(k + 1)} = \underset{f}{arg min} - μ^{T} (y - f) - ρ^{T} (z - Df) - γ^{T} (u - Ψ^{T} f) \\ + \frac{δ}{2} {∥ y - f ∥}_{2}^{2} + \frac{θ}{2} {∥ z - Df ∥}_{2}^{2} + \frac{η}{2} {∥ u - Ψ^{T} f ∥}_{2}^{2} \end{matrix}$ (11) $y^{(k + 1)} = \underset{y}{arg min} \frac{1}{2} {∥ y - g ∥}_{2}^{2} - μ^{T} (y - f) + \frac{δ}{2} {∥ y - f ∥}_{2}^{2}$ (12) $z^{(k + 1)} = \underset{z}{arg min} λ_{1} {∥ z ∥}_{1} - ρ^{T} (z - Df) + \frac{θ}{2} {∥ z - Df ∥}_{2}^{2}$ (13) $u^{(k + 1)} = \underset{u}{arg min} λ_{2} {∥ u ∥}_{1} - γ^{T} (u - Ψ^{T} f) + \frac{η}{2} {∥ u - Ψ^{T} f ∥}_{2}^{2}$ (14)

And update Lagrange multipliers: $μ^{(k + 1)} = μ^{(k)} - δ (y^{(k + 1)} - f^{(k + 1)})$ (15) $ρ^{(k + 1)} = ρ^{(k)} - θ (z^{(k + 1)} - {Df}^{(k + 1)})$ (16) $γ^{(k + 1)} = γ^{(k)} - η (u^{(k + 1)} - Ψ^{T} f^{(k + 1)})$ (17)

The sub-problem is to minimize the augmented Lagrangian function, these sub-problems are investigated one by one.

(a) f-sub-problem: Problem (11) can be solved by considering the first-order optimality condition, which can be solved by considering the following normal equation: $(δ I + θ D^{T} D + η Ψ Ψ^{T}) f^{(k + 1)} = Ψ (η u - γ) + (δ y - μ) + D^{T} (θ z - ρ)$ (18)

Since the matrix ΨΨ^T = I the matrix δI + θD^TD + ηΨΨ^T can be simplified as (δ + η) I + θD^TD . D^TD is a circulant matrix, D^TD can be diagonalized using the 2-D DFT matrix [37]. Hence, (18) has the following solution: $f^{(k + 1)} = F^{- 1} [\frac{F (Ψ (η u - γ) + (δ y - μ) + D^{T} (θ z - ρ))}{((δ + η) I + θ {| F (D) |}^{2})}]$ (19)

Where F denotes the 2D FFT transform, F^-1 denotes the 2D inverse FFT transform, and |F (D) |² denotes the magnitude square of the eigenvalues of the differential operator D.

(b) y-sub-problem: Problem (12) can be solved by considering the first-order optimality condition, which can be solved by considering the following normal equation: $(I + δ I) y^{(k + 1)} = (g + μ + δ f)$ (20) $y^{(k + 1)} = \frac{1}{1 + δ} (g + μ + δ f)$ (21)

In order to improve convergence rate, let g = f, then Equation (21) can be rewritten as shown Equation (22). $y^{(k + 1)} = \frac{1}{1 + δ} (f + μ + δ f) = f + \frac{μ}{1 + δ}$ (22)

(c) z-subproblem: Problem (13) is known as the z-sub problem, which can be solved using a shrinkage formula. The solution is given by the shrinkage formula³⁶: $z^{(k + 1)} = max {| Df + \frac{ρ}{θ} | - \frac{λ_{1}}{θ}, 0} sgn (Df + \frac{ρ}{θ})$ (23)

(d) u-sub-problem: Problem (14) is known as the u-sub problem, which can also be solved using a shrinkage formula. The solution is given by the shrinkage formula [38]: $u^{(k + 1)} = max {| Ψ^{T} f + \frac{γ}{η} | - \frac{λ_{2}}{η}, 0} sgn (Ψ^{T} f + \frac{γ}{η})$ (24)

In summary, the implementation steps of TV-L1 algorithm are given as follows:

Initialization: Input an image g; f⁽⁰⁾ = g, z⁽⁰⁾ = Dg, u⁽⁰⁾ = Ψ^Tg; The iteration number of solving the subproblems(17)–(23)is labeled by k; Regularization parameter λ₁ and λ₂; The penalty parameters δ, θ, and η; The decomposition level of wavelet function.

While stopping criteria is not met do.

Compute f subproblem using (19).

Compute y, z and u subproblem using (22),(23) and(24).

Update Lagrange multipliers by (15), (16) and (17).

End while.

Output the final correction image.

In our study, g was reconstructed using a filtered back-projection (FBP) algorithm, this TV-L1 algorithm is referred to as the FBP-TV-L1 algorithm.

2.4 FBP-GIF-TV-L1 algorithm

FBP-GIF-TV-L1 algorithm is the FBP-TV-L1 algorithm with GIF step. The FBP-GIF-TV-L1 algorithm has two main steps: GIF step and TV-L1 step.

The FBP-GIF-TV-L1 algorithm involves:

Initialization: Construct an image $f^{FBP}$ by applying filtered back-projection (FBP) algorithm; The input image f^input = f^FBP; The guidance image f^guide = f^FBP; The regularization parameter ɛ of GIF step, and the radius r of the square window of GIF step; The iteration number k of solving the subproblems (18)–(24) of TV-L1 step; Regularization parameter λ₁ and λ₂ of TV-L1 step; The penalty parameters δ, θ, and η of TV-L1 step; The decomposition level of wavelet function of TV-L1 step.

GIF step:

f^{out} = f_{boxfilter} (a) \cdot * f^{guide} + f_{boxfilter} (b)

(25)

With a and b defined as $a = cov (f^{guide} f^{input}) \cdot / (var (f^{guide}) + ɛ)$ (26) $b = f_{boxfilter} (f^{input}) - a \cdot * f_{boxfilter} (f^{guide})$ (27)

Where f_boxfilter denotes a box filter of square window. The regularization parameter ɛ in Equation (26) specify which edges or high-variance patches should be preserved.

TV-L1 step: Calculate the correction image f using TV-L1 algorithm.

min_{f} \frac{1}{2} {∥ f - f^{out} ∥}_{2}^{2} + λ_{1} {∥ f ∥}_{TV} + λ_{2} {∥ Ψ^{T} f ∥}_{1}

(28)

FBP-GIF-TV-L1 algorithm is performed by two components: GIF step, and TV-L1 step. In our study, we consider the special case f^guide = f^input = f^FBP for GIF step, then we get the output image f^out. GIF step employs one regularization parameter ɛ. TV-L1 step employs two regularization parameters: λ₁ and λ₂. The regularization parameters ɛ, λ₁, and λ₂ play an important role in ring artifacts correction. We understand that more fine tuning of the parameters may lead to better results, however, it is difficult to choose the optimal parameters. As far as we know, most parameters in regularization algorithms are chosen empirically. Therefore, in our study, we also select these parameters in an empirical fashion.

3 Numerical results

The FBP-GIF-TV-L1 algorithm can be used to correct the ring artifacts for IL-PC-CT. We compare the algorithm with four closely related methods: FBP, FBP-GIF, FBP-TV-L1, and FBP-wavelet-FFT filter. The abdomen phantom and real liver sample were used to evaluate our algorithm. All the algorithms were implemented in C++ language and matlab language on a PC (4.0 GB memory, 3.4 GHz CPU). The loss of the energy can be expressed by the energy of the difference of the reference image (t_i,j) and the corrected (r_i,j) images, relative to the reference image, resulting in the relative mean square error(MSEr). The smaller the value of MSEr, the better the corrected reconstruction results. ${MSE}_{r} = \frac{\sum_{i = 1}^{I} \sum_{j = 1}^{J} (t_{i, j} - r_{i, j})^{2}}{\sum_{i = 1}^{I} \sum_{j = 1}^{J} (t_{i, j})^{2}}$ (29)

We use the structural similarity (SSIM) index [39] as a quantitative measure for the image similarity. The computed local similarity index is defined on windows x and y: $ssim (x, y) = \frac{(2 μ_{x} μ_{y} + c_{1}) (2 σ_{xy} + c_{2})}{(μ_{x}^{2} + μ_{y}^{2} + c_{1}) (σ_{x}^{2} + σ_{y}^{2} + c_{2})}$ (30) where μ_x and μ_y are the averages of x and y; $σ_{x}^{2}$ , $σ_{y}^{2}$ are the variances; σ_xy is the covariance of x, y; and c₁, c₂ are two variables to stabilize the division with weak denominator. The overall SSIM is the mean of local similarity indices: $SSIM (X, Y) = \frac{1}{N} \sum_{i = 1}^{N} ssim (x_{i}, y_{i})$ (31) where X is an image, Y is another image; x_i, y_i are the corresponding windows indexed by i; and N is the number of windows. Here, the size of the windows is 8 × 8. The larger value of SSIM, the more similarity between two images.

3.1 The ring artifacts correction method FBP-GIF

3.1.1 Simulation data

As shown in Fig. 3 (a), the abdomen phantom without ring artifacts is used as a reference image to measure the performance of the FBP-GIF algorithm. The image size is 512×512 pixels, and 512 projections were used for these constructions. Two regions of interest (ROI) (Fig. 3(a)) containing both smooth components and details are magnified. In the simulation test, lines with variable intensity were added to the sonogram. As shown in Fig. 3(b), image reconstructed by FBP was corrupted by ring artifacts after adding variable intensity lines in the sonogram. As shown in Fig. 3(c)-(f), GIF step suppressed the ring artifacts when the parameter ɛ increases, however, some minor residual ring artifacts may still remain. Especially, as shown in Fig. 3 (e) and 3(f), the ROI containing details were over smoothed and result in a blurred image using the bigger parameter ɛ. A good image-filtering method smooths an image by suppressing noise while preserving important image details. Taking this factor into account: For simulated phantom, the parameter ɛ = 0.00935 was chosen.

Fig.3

(a) reference image. (b) FBP. (c) FBP-GIF (ɛ = 0.000935). (d) FBP-GIF (ɛ = 0.00935). (e) FBP-GIF (ɛ = 0.0935). (f) FBP-GIF (ɛ = 0.935).

To compare the results in a quantitative way, we have computed the MSEr and SSIM, which are widely used to measure image qualities. From Table 1, we see that the correction reconstruction image obtained by FBP-GIF (ɛ = 0.00935) is of the best quality.

Table 1

MSEr and SSIM (%) comparison

	FBP	FBP-GIF	FBP-GIF	FBP-GIF	FBP-GIF
	FBP	(ɛ = 0.000935)	(ɛ = 0.00935)	(ɛ = 0.0935)	(ɛ = 0.935)
MSE_r	0.0054	0.0056	0.0048	0.0078	0.0137
SSIM	96.68	97.08	97.56	94.51	93.27

3.1.2 Real data

As shown in Fig. 4 (a), the image of real liver sample reconstructed by FBP was corrupted by ring artifacts. For real liver sample, we also evaluate the ring artifacts correction performance of FBP-GIF for different regularization parameter ɛ in Fig. 4(b)-4(d). The results show that GIF step can suppress the ring artifacts when the parameter ɛ increases, however, some residual ring artifacts may still remain.

Fig.4

(a) FBP. (b) FBP-GIF(ɛ = 0.00015). (c) FBP-GIF (ɛ = 0.00075). (d) FBP-GIF (ɛ = 0.00135).

In order to evaluate the correction results objectively, as shown in Figs. 5 and 6, vertical line I(Fig. 4(a)) profiles through the center of the image show that the “high variance” of reconstructed slices affected to be small by the guided image filtering using three different parameters ɛ (Fig. 4 solid arrow). As shown in Fig. 6, vertical line II(Fig. 4(a)) profiles through the center of the image show that most of ring artifacts become piece wise smooth domain by the guided image filtering using the parameter ɛ = 0.00135 (Fig. 6 solid arrow). Namely, guided image filtering can suppress the ring artifacts to some extent. However, the parameter ɛ is not the bigger the better. Figure 7(a)-7(b) show output image using guided image filtering by two different parameters: ɛ = 0.00135 and ɛ = 0.00435. Figure 7(c)-7(d) show the corresponding magnified regions of a small region of interest. As shown in Fig. 7(c)-7(d), a larger ɛ = 0.00435 cannot preserve the liver section details (solid arrow). For real data, the parameter ɛ = 0.00135 was chosen.

Fig.5

line I (Fig. 3(a)) profile through the center of the reconstructed slices in Fig. 3.

Fig.6

line II(Fig. 3(a)) profile through the center of the reconstructed slices in Fig. 3.

Fig.7

(a) FBP-GIF(ɛ = 0.00135). (b) FBP-GIF(ɛ = 0.00435) and the corresponding magnified regions. (c) FBP-GIF(ɛ = 0.00135). (d) FBP-GIF(ɛ = 0.00435).

3.2 Comparison of ring artifact correction

3.2.1 Simulation data

The main purpose of this section is to evaluate the performance of the FBP-GIF-TV-L1, which was compared with four closely related methods: FBP, FBP-GIF, FBP-wavelet-FFT filter, and FBP- TV-L1. The parameters were set to λ₁ = 2e - 4 and λ₂ = 7e - 4 for FBP-GIF-TV-L1 algorithm. The maximum iteration number of solving the FBP-GIF-TV-L1 sub-problems performed k = 40. For wavelet function, we use Daubechie wavelet with decomposition level 5. The penalty parameters were set empirically to δ = 0.01, θ = 1.1, and η = 0.001. The parameter ɛ for guided image filtering was empirically set to ɛ = 0.00935, the box-filter-window parameter to r = 4. Wavelet-FFT filter is a fast, powerful and stable filter based on combined wavelet and FFT transforms was presented to preprocess sonogram. Wavelet-FFT filter employs three parameters for the filtering process¹⁴: the highest decomposition level L, the wavelet type and the damping factor σ. In the first step, a tight condensation of the stripe information by combining the wavelet and FFT transforms. Then it performs a damping of the relevant coefficients by using a Gaussian function. When the sinogram preprocessing is complete, the output sinogram can be used to reconstruct image by the FBP algorithm. The reconstruction algorithm based on sinogram preprocessing is referenced as FBP-wavelet-FFT filter.

As for FBP-wavelet-FFT filter, the impact of L is illustrated in Fig. 8. When we fix the damping factor σ = 7_, Fig. 8 show that the image quality is similar to each other when we choose different the highest decomposition level L from 2, 4, 6. In Fig. 9, we choose different damping factor σ from 1, 3, 5, 7, 9 when we fix the highest decomposition level L = 4. Figure 9 shows that σ = 7 should be chosen for the best visualization.

Fig.8

FBP-wavelet-FFT-filter with differernt decomposition level L.

Fig.9

FBP-wavelet-FFT-filter with differernt damping factor.

For convenience of comparison, we present the images from Fig. 10. Two regions of interest A and B are magnified in Figs. 10 and 11. The results in Fig. 10(c) and Fig. 11(c) show that image reconstructed by FBP was corrupted by ring artifacts. The result from Fig. 10(d) has been obtained with the Daubechie wavelet L = 4 and σ= 7. The results in Fig. 10(d) and Fig. 11(d) show that most ring artifacts are reduced from the reconstructed slice. However, Fig. 11(d) also shows that image reconstructed by FBP-wavelet-FFT filter still contains minimal distortion at the image center. FBP-GIF algorithm can suppress the ring artifacts, however, some ring artifacts may still remain in Fig. 10(e) and Fig. 11(e). Figure 10(b) and Fig. 11(b) show thatFBP-GIF-TV-L1 algorithm can removed the ring artifacts completely and preserve the detail.

Fig.10

(a) reference image. (b) FBP-GIF-TV-L1. (c) FBP. (d) FBP-wavelet-FFT-filter. (e) FBP-GIF.

Fig.11

Partial magnification: (a) reference image. (b) FBP-GIF-TV-L1. (c) FBP. (d) FBP-wavelet-FFT-filter. (e) FBP-GIF.

As shown in Fig. 12, vertical line I (Fig. 10(a)) profiles through the center of the image show that the gray level value of curve obtained by FBP-GIF-TV-L1 algorithm is the closest to the gray level value of reference image among the four curves. Table 2 summarizes the MSEr and SSIM values between images reconstructed using different algorithm and reference image. It is worth noticing that image reconstructed by FBP-GIF-TV-L1 better than other algorithms.

Fig.12

line I(Fig. 7(a)) profile through the center of the reconstructed slices in Fig. 7.

Table 2

MSEr and SSIM (%) comparison

	FBP	FBP-GIF-TV-L1	FBP-GIF	FBP-wavelet-
		(ɛ = 0.00935, k = 40)	(ɛ = 0.00935)	FFT-filter
MSE_r	0.0054	0.0037	0.0048	0.0184
SSIM	96.68	97.13	97.06	94.59

In order to demonstrate the advantage of the GIF step, we compare FBP-TV-L1 with or without GIF step. We also discuss the impact of the maximum iteration number k of solving the TV-L1 sub-problems for FBP-TV-L1. The image results are shown in Fig. 13. We evaluate the ring artifacts correction performance of FBP-TV-L1 for different iteration number k in Fig. 13(a) and 13(c). Particularly, we choose k from 40, 480. As shown in Fig. 13(d) and 13(f), FBP-GIF-TV-L1 (k = 40) and FBP-TV-L1 (k = 480) outperforms FBP-TV-L1 (k = 40), in terms of suppressing ring artifacts. However, the results show that some ring artifacts may still remain using FBP-TV-L1 algorithm. Figure 13(b) and 13(e) suggest that FBP-GIF-TV-L1 can remove the ring artifacts completely.

Fig.13

(a) FBP-TV-L1 (k = 40). (b) FBP-GIF-TV-L1(k = 40). (c)FBP-TV-L1(k = 480) and the corresponding magnified regions. (d) FBP-TV-L1(k = 40). (e) FBP-GIF-TV-L1(k = 40). (f) FBP-TV-L1(k = 40).

In order to evaluate the correction reconstruction results objectively, SSIM and runtime are adopted to evaluate the algorithm, just as in Table 3. When the SSIM value between reference image and image corrected using FBP-GIF-TV-L1 is larger compared with FBP-TV-L1 algorithm, Table 3 shows that FBP-GIF-TV-L1 algorithm provides considerable time benefits over the FBP-TV-L1 algorithm (125.125s using FBP-GIF-TV-L1 algorithm compared with 1076.453s using FBP-TV-L1 algorithm (k = 480)). Moreover, the SSIM also demonstrate that the FBP-GIF-TV-L1 can produce higher accuracy results.

Table 3

Runtime (sec) and SSIM (%) comparison

	FBP-TV-L1	FBP-GIF-TV-L1	FBP-TV-L1
	(k = 40)	(k = 40)	(k = 480)
TIME	124.508	125.125	1076.453
SSIM	96.75	97.13	97.08

3.2.2 Real data

First, to show the practical applicability of our approach, we test our proposed algorithm on the liver sample. Figure 16 shows the results obtained using four different algorithms: FBP, FBP-GIF-TV-L1, FBP-wavelet-FFT filter, and FBP-GIF. The parameters λ_i (i = 1, 2) of the regularization parameters were set empirically to λ₁ = 2e - 8 and λ₂ = 7e - 3 for FBP-GIF-TV-L1 algorithm. The maximum iteration number of solving the sub-problems (18)–(24) performed k = 10. For wavelet function, we use Daubechie wavelet with decomposition level 5. The parameters δ, θ, and η of the penalty parameters were set empirically to δ = 0.01, θ = 1.1, and η = 0.001. As for FBP-wavelet-FFT filter, the impact of L is illustrated in Fig. 14. When we fix the damping factor σ = 2.4, Fig. 14 show that the image quality is similar to each other when we choose different the highest decomposition level L from 4, 5, 6. In Fig. 15, we choose different damping factor σ from 0.4, 2.4, 4.4 when we fix the highest decomposition level L = 5. Figure 15 shows that σ = 2.4 should be chosen for the best visualization. The parameter ɛ for guided image filtering was empirically set to ɛ = 0.00135, the box-filter-window parameter was set to r = 4.

Fig.14

FBP-wavelet-FFT-filter with differernt decomposition level L.

Fig.15

FBP-wavelet-FFT-filter with differernt damping factor.

Fig.16

(a) FBP. (b) FBP-GIF-TV-L1. (c) FBP-wavelet-FFT filter. (d) FBP-GIF.

Figures 16(a) and 17(a) show that the ring artifacts generated by the FBP algorithm. Figures 16(d) and 17(d) show that FBP-GIF algorithm can suppress the ring artifacts to some extent, however, some residual ring artifacts may still remain. The result from Fig. 16(c) has been obtained with Daubechie wavelet L = 5 and σ = 2.4. Figures 16(c) and 17(c) show that the ring artifacts are suppressed well by the FBP-wavelet-FFT filter algorithm, with minimal distortion at the image center (Fig. 17(c) solid arrow). The experiment result in Figs. 16(b) and 17(b) suggests that FBP-GIF-TV-L1 algorithm can removed the ring artifacts completely and preserve the soft tissue detail (solid arrow). Besides, additional artifacts are less likely to appear thanks to the regularization. Figure 17(b) shows that the image reconstructed by FBP-GIF-TV-L1 method is visually more satisfying, and succeeds in removing almost all visible rings to an appreciable level.

Fig.17

Partial magnification. (a) FBP. (b) FBP-GIF-TV-L1. (c) FBP-wavelet-FFT filter. (d) FBP-GIF.

Next, for real liver sample, we also evaluate the ring artifacts correction performance of FBP-TV-L1 for different iteration number k. The images were shown in Fig. 18 (a) and 18(c). Figure 18(a) show that FBP-TV-L1 (k = 10) algorithm can suppress the ring artifacts to some extent, however, some ring artifacts may still remain. In order to demonstrate the effective and computationally efficient of the GIF step, we compare FBP-TV-L1 with or without GIF step. Figure 18(b)-(c) suggest that FBP-GIF-TV-L1 and FBP-TV-L1 (k = 40) algorithm can removed the ring artifacts completely and preserve the soft tissue detail (solid arrow).

Fig.18

(a) FBP-TV-L1(k = 10); (b) FBP-GIF-TV-L1(k = 10);(c) FBP-TV-L1(k = 40). andthe corresponding magnified regions (d)magnified ROI of FBP-TV-L1(k = 10). (e)magnified ROI of FBP-GIF-TV-L1(k = 10); (f)magnified ROI of FBP-TV-L1(k = 40).

We also discuss the impact of the maximum iteration number k of solving the TV-L1 sub-problems for FBP-TV-L1. In order to make a comparison study, one regions of interest were magnified. It is clear from Fig. 18(d) and 18(f) that FBP-TV-L1 (k = 40) outperforms FBP-TV-L1 (k = 10) in terms of suppressing ring artifacts. Figure 18 (e) and 18 (f) show that FBP-GIF-TV-L1 and FBP-TV-L1 (k = 40) can remove the ring artifacts completely and preserve the detail. For the convenience of discussion, FBP-GIF-TV-L1 (k = 10) image is used as a reference image. When the image similarity value SSIM = 99.7% between FBP-GIF-TV-L1 (k = 10) image and FBP-TV-L1 (k = 40) image, Table 4 shows that the computational time of FBP-TV-L1 (k = 40) increases compared with FBP-GIF-TV-L1.

Table 4

Runtime (sec) and SSIM (%) comparison

	FBP-TV-L1	FBP-GIF-TV-L1	FBP-TV-L1
	(k = 10)	(ɛ = 0.00135, k = 10)	(k = 40)
TIME	333.89	336.02	647.19
SSIM	97.80	100	99.70

4 Conclusion

Our proposed a new FBP-GIF-TV-L1 algorithm to correct the ring artifacts and preserve fine structures for IL-PC-CT. The new correction method has two main steps: the GIF step and the TV-L1 step. The GIF step can suppress the input image f^FBP corrupted by ring artifacts to some extent. Then the TV-L1 Model further effectively correct ring artifacts. We apply ADMM to solve the TV-L1 Model.

To evaluate the proposed algorithm practicality and efficiency, the simulation data experiments and real synchrotron data experiments are conducted. For comparison, the FBP-wavelet-FFT filter is used to correct ring artifacts on the same data set. Based on the simulations and real synchrotron data experiments, we come up with the following conclusions. (1) Both FBP-GIF-TV-L1 and FBP-TV-L1 algorithm are effective in correcting the wide and faint ring artifacts. (2) The FBP-GIF-TV-L1 algorithm is an accelerated version of the FBP-TV-L1 algorithm in removing almost all visible rings to an appreciable level. (3) FBP-GIF-TV-L1 algorithm can improve the image quality over the FBP-wavelet-FFT filter in terms of preserving characteristics of the soft tissue and reducing the ring artifacts.

However, the real liver sample is corrupted by wide and faint ring artifacts, we have only studied that the FBP-GIF-TV-L1 algorithm is effective for the wide and faint ring artifacts. Further studies are currently under way to correct strong ring artifacts.

Footnotes

Acknowledgments

This work was supported by National Natural Science Foundation of China (61671004, 61271012, 81371549, 81671683, and 11501415); the Natural Science Foundation of Tianjin City in China (16JCYBJC28600); the WBE Liver Fibrosis Foundation of China (Grant No.CFHPC20131033); the Instrument Developing Project of the Chinese Academy of Sciences (YZ201410); the Foundation of Tianjin university of technology and education (KJ11-22; KJ17-36; J10011060321); The IHEP-CAS Scientific Research Foundation (2013IHEPYJRC801); The Open Project of Key laboratory of Opto-electronic Information Technology, Ministry of Education (Tianjin University) (Grant No. 2017KFKT004); The Foundation of Tianjin Municipal Education Commission (JWK1603).

References

Boin

and Haibel

, Compensation of ring artefacts in synchrotron tomographic images, Optics Express 14 (2006), 12071–12075.

Gao

, Luo

, Yin

, et al., A micro-tomography method based on X-ray diffraction enhanced imaging for the visualization of micro-organs and soft tissues, Computerized Medical Imaging & Graphics 30 (2006), 339–347.

Mayo

, Davis

, Gureyev

T.P.

, et al., X-ray phase-contrast microscopy and microtomography, Optics Express 11 (2003), 2289–2302.

Chen

R.C.

, Rigon

and Longo

, Quantitative 3D refractive index decrement reconstruction using single-distance phase-contrast tomography data, Journal of Physics D Applied Physics 44 (2011), 495401:1–495401:9.

Liu

, Cai

and Ning

, Evaluation of differential phase contrast cone beam CT imaging system, Journal of X-ray Science and Technology 25 (2017), 357–372.

Danon

, Block

R.C.

and Slovacek

R.E.

, Design and construction of a thermal neutron target for the RPI linac, Nuclear Instruments & Methods in Physics Research 352 (1995), 596–604.

Dilmanian

, Zhong

, Ren

, Wu

, et al., Computed tomography of x-ray index of refraction using the diffraction enhanced imaging method, Physics in Medicine & Biology 45 (2000), 933–946.

Wang

, Huang

, Chen

, et al., Low-dose multiple-information retrieval algorithm for X-ray grating-based imaging, Nuclear Inst & Methods in Physics Research A 635 (2001), 103–107.

Pfeiffer

, Kottler

, Bunk

and David

, Hard x-ray phase tomography with low-brilliance sources, Physical Review Letters 98 (1081), 108105:1–108105:4.

10.

Meiser

, Willner

, Schroter

, et al., Increasing the field of view in grating based X-ray phase contrast imaging using stitched gratings, Journal of X-ray Science and Technology 24 (2016), 379–388.

11.

Raven

, Snigirev

, Snigireva

and Spanne

, Phase-contrast Microtomography with Coherent High-energy Synchrotron X-rays, Applied Physics Letters 69 (1996), 1826–1828.

12.

, Hu

and Yang

, Image reconstruction algorithm for in-line phase contrast imaging computed tomography with an improved anisotropic diffusion model, Journal of X-ray science and technology 23 (2015), 311–320.

13.

Herman

G.T.

and Brouw

W.N.

, Book-review - image reconstruction from projections - the fundamentals of computerized tomography, 29 (1980), 64–68.

14.

Davis

G.R.

and Elliott

J.C.

, X-ray microtomography scanner using time-delay integration for elimination of ring artefacts in the reconstructed image, Nuclear Instruments & Methods in Physics Research 394 (1997), 157–162.

15.

Görner

, Hentschel

M.P.

, Müller

B.R.

, et al., BAMline: The first hard X-ray beamline at BESSY II, Nuclear Instruments & Methods in Physics Research s 467-468 (2001), 703–706.

16.

Münch

, Trtik

, Marone

and Stampanoni

, Stripe and ring artifact removal with combined wavelet–Fourier filtering, Optics Express 17 (2009), 8567–8591.

17.

Huo

Q.R.

, Li

J.W.

and Lu

, Removing ring artifacts in CBCT images via L0 smoothing, Imaging Systems and Technology 26 (2016), 284–294.

18.

Sijbers

and Postnov

, Reduction of ring artefacts in high resolution micro-CT reconstructions, Physics in Medicine & Biology 49 (2004), 247–253.

19.

Prell

, Kyriakou

W.A.

and Kalender

, Comparison of ring artifact correction methods for flat-detector CT, Physics in Medicine & Biology 54 (2009), 3881–3895.

20.

Yuan

, Qi

. Chen

, et al, CT image ring artifact reduction based on prior image, Med Phys 43 (2016), 3388.

21.

Pierre

and Alessandro

, Ring artifacts correction in compressed sensing tomographic reconstruction, Journal of Synchrotron Radiation 22 (2015), 1268–1295.

22.

Raven

, “Numerical removal of ring artifacts in microtomography, Review of Scientific Instruments 69 (1998), 2978–2980.

23.

, Sun

and Tang

, Guided image filtering, IEEE Transactions on Pattern Analysis and Machine Intelligence 35 (2013), 1397–1409.

24.

, Qu

and Liu

, Simultaneous algebraic reconstruction technique based on guided image filtering, Optics Express 24 (2016), 15897–15911.

25.

Tibshirani

R.J.

, Regression shrinkage and selection via the LASSO. J R Stat Soc B, Journal of the Royal Statistical Society 58 (1996), 267–288.

26.

Shi

, et al., Enhanced spatial resolution in fluorescence molecular tomography using restarted l1-regularized nonlinear conjugate gradient algorithm, Journal of Biomedical Optics 19 (2014), 046018:1–046018:11.

27.

Laroque

S.J.

, Sidky

E.Y.

and Pan

, Accurate image reconstruction from sparse data in diffraction tomography using a total variation minimization algorithm, Proc SPIE 6513 (2017), 651302–651306.

28.

Wang

and Yu

, A soft-threshold filtering approach for reconstruction from a limited number of projections, Physics in Medicine and Biology 55 (2010), 3905–3916.

29.

D.J.

, Qu

G.R.

, Hu

C.H.

, et al., Anisotropic total variation minimization approach in in-line phase-contrast tomography and its application to correction of ring artifacts, Chinese Physics B 26 (2017), 060701:1–060701:8.

30.

Gao

and Zhao

, Multilevel bioluminescence tomography based on radiative transfer equation part 2: Total variation and l1 data fidelity, Optics Express 18 (2010), 1854–1871.

31.

Liu

L.K.

, Chan

S.H.

and Nguyen

T.Q.

, Depth reconstruction from sparse samples: Representation, algorithm, and sampling, IEEE Transactions on Image Processing 24 (2015), 1983–1996.

32.

Draper

and Smith

, Applied Regression Analysis (John Wiley, 1981).

33.

Lauzier

P.T.

and Chen

G.H.

, Characterization of statistical prior image constrained compressed sensing (PICCS): II. Application to dose reduction, Medical Physics 40 (2013), 021902:1–021902:14.

34.

Jiang

, Liu

, An

, et al., Novel l 2,1-norm optimization method for fluorescence molecular tomography reconstruction,:, Biomedical Optics Express 7 (2016), 2342–2359.

35.

Correia

, et al., Wavelet-based data and solution compression for efficient image reconstruction in fluorescence diffuse optical tomography, Journal of Biomedical Optics 18 (2013), 086008:1–086008:13.

36.

Boyd

, Parikh

, Chu

, et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations & Trends ® in Machine Learning 3 (2011), 1–122.

37.

M.K.

, Chan

R.H.

and Tang

W.C.

, A Fast Algorithm for Deblurring Models with Neumann Boundary Conditions, Siam Journal on Scientific Computing 21 (1999), 851–866.

38.

, An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing, Dissertations & Theses - Gradworks (2011).

39.

Wang

, Bovik

A.C.

, Sheikh

H.R.

and Simoncelli

E.P.

, Image quality assessment: From error visibility to structural similarity, IEEE Trans. Image Process 13 (2004), 600–612.