First-order primal-dual algorithm for sparse-view neutron computed tomography-based three-dimensional image reconstruction

ACCELERATOR, RAY TECHNOLOGY AND APPLICATIONS

First-order primal-dual algorithm for sparse-view neutron computed tomography-based three-dimensional image reconstruction

Yang Liu ，

Teng-Fei Zhu，

Zhi Luo，

Xiao-Ping Ouyang

Nuclear Science and Techniques

Vol.34, No.8

Article number 118

Published in print Aug 2023

Available online 21 Aug 2023

DOI：10.1007/s41365-023-01258-4

54400

Neutron computed tomography (NCT) is widely used as a noninvasive measurement technique in nuclear engineering, thermal hydraulics, and cultural heritage. The neutron source intensity of NCT is usually low and the scan time is long, resulting in a projection image containing severe noise. To reduce the scanning time and increase the image reconstruction quality, an effective reconstruction algorithm must be selected. In CT image reconstruction, the reconstruction algorithms used can be divided into three categories: analytical algorithms, iterative algorithms, and deep learning. Because the analytical algorithm requires complete projection data, it is not suitable for reconstruction in harsh environments, such as strong radiation, high temperature, and high pressure. Deep learning requires large amounts of data and complex models, which cannot be easily deployed, as well as has a high computational complexity and poor interpretability. Therefore, this paper proposes the OS-SART-PDTV iterative algorithm, which uses the ordered subset simultaneous algebraic reconstruction technique (OS-SART) algorithm to reconstruct the image and the first-order primal-dual algorithm to solve the total variation (PDTV), for sparse-view NCT three-dimensional reconstruction. The novel algorithm was compared with other algorithms (FBP, OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV) by simulating the experimental data and actual neutron projection experiments. The reconstruction results demonstrate that the proposed algorithm outperforms the FBP, OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms in terms of preserving edge structure, denoising, and suppressing artifacts.

NCTFirst-order primal-dual algorithmOS-SARTTotal variationSparse-view

Introduction

Neutron tomography is a noninvasive measurement method that differs from X-ray imaging, and muon tomography has been applied in various fields with a high accuracy and reliability [1,2]. The strong transmission of neutrons in metallic elements makes neutron tomography extremely useful for examining hydrogen in heavy elements under harsh experimental conditions such as strong radiation, high-temperature, and high-pressure two-phase flow systems. In addition, neutrons are used in boron neutron therapy owing to their ability to detect and reconstruct the three-dimensional (3D) distribution of boron concentration [3,4]. Neutron tomography systems consist of two parts: a neutron photographic system and a computed tomography (CT) system capable of measuring and displaying the 3D void rate in the boiling flow in a heated rod beam. They have also been applied to simulate the cores of advanced nuclear [5-8]. Asano and Takenaka [9] used neutron tomography to determine the void-rate distribution in two-phase air and water flows. Kureta et al. visualized the toroidal streaming in a heated fuel rod bundle using 3D neutron tomography [10]. Tremsin et al. used neutron imaging to assess the attachment, deposition, and structural integrity of fuel core blocks [11]. Andersson et al. developed a portable fast NCT to obtain the void rate distribution in fuel rod bundles in a boiling water reactor [12]. Although portable NCT systems can perform measurements at any altitude, their neutron generators exhibit extremely low neutron yields. Therefore, promoting the rapid development of neutron tomography for improving the efficiency and safety of nuclear fuels is an important research direction.

The most crucial part of a neutron tomography system is the 3D reconstruction algorithm, which is directly related to the image quality and measurement reliability. Because of the long scanning time of neutron CT (ranging from 30 min to several hours) and the complexity of practical application environments, it is difficult to obtain complete projection data. Consequently, reconstruction algorithms that can efficiently utilize sparse-view projection data are required to obtain accurate measurement results. Furthermore, the reconstruction of sparse-view projection data is helpful for investigations in harsh environments such as strong radiation. Currently, algorithms are used for the reconstruction of sparse views and highly under-sampled projection data. Therefore, this study was focused on a 3D reconstruction algorithm for sparse-view NCT.

Sparse-view reconstruction is a highly complex mathematical problem and is a significant research direction in the image reconstruction field. Because prior information can be embedded in the iterative reconstruction algorithm, it is considered the sparse-view reconstruction algorithm most likely to replace the FBP algorithm. With the development of deep learning, it has been applied to the field of tomographic reconstruction. Many researchers have worked on sparse-view tomography and achieved significant results. As proposed in [13], the DEAR model, which adds a priori information and data in the image domain to the compressed sensing-based variation model, can not only eliminate the artifacts of the reconstructed image but also improve the quality of the reconstructed image. A multiple adversarial learning angiography image reconstruction framework has been proposed in the literature [14], which solves the challenge of low-intensity aortic reconstruction by introducing dual-correlation constrained adversarial learning; its application in clinical data shows its feasibility and effectiveness. However, because of practical limitations, this study focused on neutron CT 3D reconstruction, and it was not possible to obtain a large number of neutron projection images for supervised training. Therefore, this study focuses on an iterative algorithm for neutron CT 3D reconstruction. Currently, mainstream iterative reconstruction algorithms include ART [15], SIRT [16], and SART [17]. Although all these iterative algorithms can reconstruct high-quality images, when the projection data are extremely sparse and no a priori information is introduced, the reconstructed images show significant artifacts. Donoho [18] proposed the compressive sensing (CS) theory, which has promoted the rapid development of sparse-view image reconstruction. According to the mathematical implications of the CS theory, an image may be reconstructed exactly if it is sparse or can be sparsely represented by a spatial transformation. The evolution of the CS theory has facilitated the application of total variational algorithms in sparse-view image reconstruction. Rudin proposed the ROF denoising model, which reduces noise while preserving the image edge and detailed structural information using the total variance as a constraint for image denoising [19]. Several TV-based techniques and variations to improve the performance of the algorithm have been put forward, including directional TV, weighted TV [20], edge-preserving TV [22], and weight TV[23]. The total variation model is a typical regularization model. Total variance minimization is particularly important in CT imaging because of its ability to resolve acute discontinuities. This is critical for many imaging problems, because the sample edges contain critical structural information about the object. However, because of the non-smooth character of the total variance, it is difficult to achieve total variance regularization minimization. Furthermore, one drawback of these models is that they assume that the image is piecewise constant, which destroys significant information in the image. The total-variational image-denoising problem is a class of convex optimization problems. To solve convex optimization problems, the second-order method exhibits good convergence and requires only a small number of iterations. However, each iteration of the second-order method is extremely complicated, making it difficult to apply to large-scale problems. In contrast, the first-order method involves only function values and gradient information, and the process of each iteration is relatively simple. In the case of gradient magnitude sparsity, the total variance of the intermediate images is minimized, resulting in high-quality reconstructed images. Therefore, this study employs a first-order method to address the total-variance image noise reduction problem, and proposes an algorithm applicable to the 3D reconstruction of NCT in sparse views.

The goal of this study was to find an efficient 3D reconstruction algorithm for sparse-view NCT. To address artifacts and noise in sparse-view projection-reconstructed images, we propose the ordered subset simultaneous algebraic reconstruction technique (OS-SART-PDTV) algorithm, which uses the ordered subset simultaneous algebraic reconstruction technique (OS-SART) to reconstruct the image and the first-order primal-dual algorithm to solve the total variation (PDTV). This novel algorithm achieves a blurred version of a piecewise constant object through phenomenological modeling, and reconstructs high-quality 3D images by minimizing the total variance of the intermediate images through gradient-amplitude sparsity, which allows the reconstructed image interiors to change quickly and smoothly. Next, we present the proposed method, including the proposed algorithm and NCT fundamentals, in Section 2. In Section 3, the simulation experiments and a real neutron projection experiment are discussed. Section 4 discusses the research process and presents conclusions. Finally, a short conclusion is presented in Section 5.

Principles and Algorithms of NCT

2.1

NCT system

This section describes the fundamentals of CT image reconstruction. In general, CT image reconstruction generates a 3D volume from a 2D collection of projected images using mathematical algorithms. This requires the development of a measurement model using mathematical methods to relate the measured data to the desired physical properties. Image reconstruction consists of two processes: forward and back projection. Forward projection enables the calculation of model measurements corresponding to the actual measurements that match the physical properties. Model measurements corresponding to the actual measurements were calculated using forward projection. Forward projection refers to the mapping of an objective onto the projection domain. Back projection is the opposite of forward projection, which determines the physical properties of a measurement. Therefore, these two operations largely determine the image reconstruction precision. The change in forward projection depends on the radiation beam geometry. Based on the selected system, source of radiation, and available data acquisition systems, the appropriate geometry of the beam was selected so that reliable measurements could be obtained. The NCT system primarily includes a collimator and shield, rotating sample stage, and neutron imaging system. The neutron imaging system included a 2D neutron image detector, conversion screen, and 3D reconstruction software. The algorithm used for 3D reconstruction directly determines the imaging quality and measurement reliability of the system, and is the core of the system. A schematic of the NCT scanning system is shown in Figure 1.

Fig. 1

(Color online) Schematic of the NCT system

In NCT systems, the mass of the neutron source and noise generated by the electronics can cause the projection data to contain noise. Therefore, Poisson noise and Gaussian noise are added to the numerical projection data reconstruction to validate the denoising ability of the algorithm while testing its reliability [24]. Therefore, these two noise components must be considered during NCT. The equation is as follows: $\hat{I} = Poisson (I) + Gaussian (m_{ic} + δ_{ic}^{2}),$ (1) where $\hat{I}$ denotes noisy transmission data, I represents the mean photon number, m_ic indicates the average value of the electronic noise, and $δ_{ic}^{2}$ is the variance of the electronic noise. System calibration typically has a value from m_ic to zero. The variance $δ_{ic}^{2}$ was estimated via dark current measurements [25].

2.2

Principle of CT

2.2.1

CT projection acquisition

In NCT, the geometry of the neutron beam is primarily parallel. Among the mathematical concepts, the parallel geometric beam is the simplest. An object can be represented by a binary function $R (x, y)$ , which represents the intensity function on the x-y plane. The projection represents the decay integral of the neutron beam as it passes through the object from the neutron source, which is considered to be the sum of the line integral and ray. The expression of the projection is given by the equation $h (s, β) = \int_{L} R (x, y) d s,$ (2) where $L$ represents a neutron beam passing straight through an object $R (x, y)$ , $s$ represents the detector coordinates in the projection domain, and $β$ denotes the counterclockwise rotation angle of the x-axis in the projection domain. According to the Delta functions in Eq. (2) can be rewritten as follows: $h (s, β) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} R (x, y) δ (x \cos β + y \sin β - s) d x d y .$ (3)

Eq. (3) represents the Radon transform of the object $R (x, y)$ , which transforms the space data into the frequency domain. The projection consisted of a group of integrals of the neutron attenuation coefficients of the sample measured at the detector for every given angle.

2.2.2

Central slice theorem

The Central slice theorem establishes a connection between the 2D projection of an object and its 1D projection in the Fourier domain. The 1D Fourier transform function $K (σ)$ of the projection function $h (s, β)$ of the density function $R (x, y)$ in a certain direction is the value of the 2D Fourier transform of the density function $R (x, y)$ on the $x - y$ plane along the same direction on a straight line past the origin. Eq. (4) represents the defined 1D Fourier variation of the objective function in polar coordinates. $K (σ, β) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} R (x, y) e^{- j 2 π σ (x \cos β + y \sin β)} d x d y$ (4)

Eq. (4) represents the 2D Fourier transform of the spatial frequency $σ_{x} = σ \cos β, σ_{y} = σ \sin β$ in the Fourier space. Therefore, the following equation was derived: $R (x, y) = \int_{0}^{2 π} \int_{0}^{\infty} Q (σ_{x}, σ_{y}) e^{j 2 π σ (x \cos β + y \sin β)} d σ d β .$ (5)

Several problems arise when the central slice theorem is used for the Fourier transform reconstruction. First, the theorem produces data in the Fourier space that are inconsistent with the Cartesian space. The data were interpolated into Cartesian coordinates. However, interpolation in the Fourier space before the anti-Fourier transform can have a significant influence on the reconstruction. The second problem is the difficulty in performing targeted reconstruction and identifying fine structures within small areas.

2.3

Algorithms

2.3.1

TV algorithm

In the NCT, the measurement time can be drastically reduced by reducing the number of projection views; however, this leads to sparse sampling. Based on the CS theory, a signal can be reconstructed exactly if it is sparse or can be represented sparsely. Under ideal conditions, NCT reconstruction can be converted to solve the problem of a linear system; thus, the discrete model of CS theory may be represented. $L g = b,$ (6) where $L = ｛ l_{i, j} ｝$ stands for system matrix $b$ denotes the projection data obtained at each angle and $g$ represents the image to be reconstructed. Mathematically, solving Eq. (6) usually converts it into a least square problem: $\arg \min_{g} {‖ L g - b ‖}_{2}^{2} .$ (7)

However, in sparse-view image reconstruction, the incompleteness of the projection data leads to reconstructed images containing noise and artifacts. Therefore, various regularization models have been proposed to obtain approximate solutions, which are defined by the following equations: $\arg \min_{g} {‖ L g - b ‖}_{2}^{2} + λ Z (g),$ (8) where $\arg \min_{g} {‖ L g - b ‖}_{2}^{2}$ represents the data fidelity term, $λ Z (g)$ denotes the regular term, and $λ$ is the regularization parameter used to balance the regular and data terms. The equation that defines TV as regularization is given in Eq. (9) ${‖ g ‖}_{TV} = \sum \sqrt{{| \nabla_{x} g |}^{2} + {| \nabla_{y} g |}^{2} + {| \nabla_{z} g |}^{2}},$ (9) where $g$ denotes the 3D image, $\nabla_{i} g$ represents the gradient along $i$ direction, and $| • |$ denotes the complex modulus.

2.3.2

First-order primal-dual algorithm

Owing to the segmental smoothness of the regularization term of the variational model, the total variational model can effectively preserve image edges during image recovery and achieve improved recovery results. Therefore, the total variational model, particularly the total variational regularization model proposed by Rudin, Osher, and Fatemi, has gained popularity. According to the ROF model, the discrete total variational model of the 3D images can be formulated as follows [26]: $\min_{g} \frac{1}{2} {‖ H g - u ‖}_{2}^{2} + α C T V (g),$ (10) where ${‖ • ‖}_{2}$ stands for the Euclidean parametrization, $α > 0$ is the regularization argument, and $C T V (g)$ is the regular term of the discrete total variational model. To obtain a stable solution to the problem of minimization, $\nabla$ , which denotes the difference algorithm on the space $ℝ^{3 n^{2}}$ , is defined as shown in Eq. (11) below: ${(\nabla g)}_{e, u, q} = ({(\nabla g)}_{e, u, q}^{h}, {(\nabla g)}_{e, u, q}^{v}, {(\nabla g)}_{e, u, q}^{w})$ (11) where $\begin{array}{l} {(\nabla g)}_{e, u, q}^{h} = ｛ \begin{array}{l} \begin{array}{l} g_{e + 1, u, q} - g_{e, u, q} e < N \\ 0, e = N \end{array} & (e = 0, 1, 2, ..., N) \end{array} \\ {(\nabla g)}_{e, u, q}^{v} = ｛ \begin{array}{l} \begin{array}{l} g_{e, u + 1, q} - g_{e, u, q} u < M \\ 0, u = M \end{array} & (u = 0, 1, 2, ..., M) \end{array} \\ {(\nabla g)}_{e, u, q}^{w} = ｛ \begin{array}{l} \begin{array}{l} g_{e, u + 1, q} - g_{e, u, q} q < Z \\ 0, q = Z \end{array} & (q = 0, 1, 2, ..., Z) \end{array} \end{array}$ g_e,u,q denotes the $((N * M) * q + e + N * u)$ -th element of $g$ , whose position in the image is $(e, u, q)$ ; ${(\nabla g)}_{e, u, q}^{h}$ , ${(\nabla g)}_{e, u, q}^{v}$ , and ${(\nabla g)}_{e, u, q}^{w}$ are the different operators in the $x$ , $y$ , and $z$ directions at a pixel $((N * M) * q + e + N * u)$ . The expression of the regular terms of the total variational model is shown below. $\begin{matrix} C T V (g) = {‖ \nabla g ‖}_{1} \\ = \sum_{0 \leq e \leq N, 0 \leq u \leq M, 0 \leq q \leq Z} \sqrt{\sum {| {(\nabla g)}_{e, u, q}^{h} |}^{2} + {| {(\nabla g)}_{e, u, q}^{v} |}^{2} + {| {(\nabla g)}_{e, u, q}^{w} |}^{2}} \end{matrix}$ (12)

In this study, the above model was improved using the finite difference matrix $R$ instead of the difference operator $\nabla$ . $C T V (g) = {‖ R g ‖}_{1}$ (13)

Thus, the expression of the model investigated in the present study can be obtained as follows: $\min_{g} \frac{1}{2} {‖ H g - u ‖}_{2}^{2} + α {‖ R g ‖}_{1} .$ (14)

Although the total variational model efficiently preserves image edges, it is difficult to obtain the minimum value of the total variational model by solving Eq. (14) directly because of the non-smooth nature of the total variational model [27]. Because the first-order primal-dual algorithm is an efficient method for solving non-smooth convex optimizations of images, we employed it to solve Eq. (14).

Definition 1: Assume that the binary function $J (χ, ν) (χ \in Ν^{n}, ν \in Ν^{n})$ is convex for $χ$ , and concave for $ν$ . If, for any $(χ \in Ν^{n}, ν \in Ν^{n})$ , there exists an $(χ^{*}, ν^{*})$ such that Eq. (15) holds, $(χ^{*}, ν^{*})$ is the saddle point of the function $J (χ, ν)$ . $J (χ^{*}, ν) \leq J (χ^{*}, ν^{*}) \leq J (χ, ν^{*})$ (15)

Chambolle and Pock suggested a first-order primal-dual algorithm to solve the following model based on the properties of the saddle points [28]: $\min_{χ \subset Ν^{n}} \max_{ν \subset Ν^{n}} £ (χ, ν) \equiv ο (χ) + 〈 K χ, ν 〉 - ϑ (ν),$ (16) where $ο$ and $ϑ$ are lower semi-continuous convex functions, and $K$ denotes the line operator. The first-order primal-dual iterations were as follow[29,30]. $χ^{(d + 1)} = \arg \min_{x \subset Ν^{n}} ο (χ) + 〈 K χ, ν 〉 - \frac{1}{2 s} {‖ χ - χ^{(d)} ‖}_{2}^{2},$ (17) $\bar{χ}^{(d + 1)} = χ^{(d + 1)} + θ (χ^{(d + 1)} - χ^{(d)}),$ (18) $ν^{(d + 1)} = \arg \min_{ν \subset Ν^{n}} 〈 K \bar{χ}^{(d + 1)}, ν 〉 - ϑ (ν) - \frac{1}{2 t} {‖ ν - ν^{(d)} ‖}_{2}^{2},$ (19) where the parameter s, t > 0 is the iteration steps for the primal and dual variables, separately; θ denotes the combinational parameter, which is used to ensure the algorithm convergence.

Definition 2: Assuming that function $ο (g) : Ν^{n} \to Ν$ , we call the formula in Eq. (20) as a conjugate function of the function $ο (g)$ [31]. $ο^{*} (ξ) = \sup_{x \subset Ν^{n}} ｛ ξ^{T} g - ο (g) ｝$ (20) where $ξ$ represents the dual vector of $g$ , $Ν$ denotes the real number space, $Ν^{n}$ represents the $k$ -dimensional real space, and $T$ denotes the transpose.

According to Definition 2, the ${‖ R g ‖}_{1}$ can be expressed as $\max_{‖ ζ ‖ \leq 1} ξ^{T} R g$ . Then, Eq. (14) can be converted into a minimal–extreme problem [32] as follows: $\min_{g} \max_{{‖ ξ ‖}_{\infty} \leq 1} \frac{1}{2} {‖ H g - u ‖}_{2}^{2} + α ξ^{T} R g .$ (21)

The iterative formula in Eq. (21) can be obtained by using a first-order primal-dual algorithm. $g^{(d + 1)} = \arg \max_{g \subset Ν^{3 n^{2}}} \frac{1}{2} {‖ H g - u ‖}_{2}^{2} + α g^{T} R^{T} ξ^{(d)} + \frac{1}{2 s} {‖ g - g^{(d)} ‖}_{2}^{2}$ (22) $\bar{g}^{(d + 1)} = g^{(d + 1)} + θ (g^{(d + 1)} - g^{(d)})$ (23) $ξ^{(d + 1)} = \arg \max_{{‖ ξ ‖}_{\infty} \leq 1} α ξ^{T} R \bar{g}^{(d + 1)} - \frac{1}{2 t} {‖ ξ - ξ^{(d)} ‖}_{2}^{2}$ (24)

The minimization problem for the subproblem g can be converted into solving the solution of the linear Eq. (25). H is a $3 \times 3$ a chunking matrix. ${(H^{T} H + \frac{1}{s} I)}^{- 1} g^{(d + 1)} = (H^{T} u + \frac{1}{s} g^{(d)} - α R^{T} ξ^{(d)})$ (25) where H^T is the transpose of H and R^T is the transpose of R. $g^{(d + 1)} = {(H^{T} H + \frac{1}{s} I)}^{- 1} (H^{T} u + \frac{1}{s} g^{(d)} - α R^{T} ξ^{(d)})$ (26)

For subproblem $ξ$ , the maximum problem of $ξ$ is first transformed into a minimization problem. $\begin{matrix} ξ^{(d + 1)} = \arg \max_{{‖ ξ ‖}_{\infty} \leq 1} α ξ^{T} R \bar{g}^{(d + 1)} - \frac{1}{2 t} {‖ ξ - ξ^{(d)} ‖}_{2}^{2} \\ = \arg \min_{{‖ ξ ‖}_{\infty} \leq 1} \frac{1}{2} {‖ ξ - (ξ^{(d)} - t α R \bar{g}^{(d + 1)}) ‖}_{2}^{2} \end{matrix}$ (27)

Then, the equation Eq. (27) is the $ξ^{(d + 1)} = P_{A} (ξ^{(d)} - t α R \bar{g}^{(d + 1)})$ , which denotes the projection of $ξ^{(d)} - t α R \bar{g}^{(d + 1)}$ onto set $A$ . $A \equiv ｛ ξ \in R^{3 n^{2}} \times R^{3 n^{2}} : {‖ ξ ‖}_{\infty} \leq 1 ｝$ (28)

The defining equations for the projection operator of any vector $q$ on set $A$ are given by Eqs. (29) (table 1) $P_{A} (q) \equiv \arg \min_{ξ \in A} {‖ ξ - q ‖}_{2}^{2}, \forall q$ (29)

First-order primal-dual algorithm

Set

s > 0, t > 0, θ \in [0, 1], input, u, α, R

Initialize

ξ^{(0)}, g^{(0)}

Repeat

k \leftarrow k + 1

g^{(d + 1)} = {(H^{T} H + \frac{1}{s} I)}^{- 1} (H^{T} u + \frac{1}{s} g^{(d)} - α R^{T} ξ^{(d)})

\bar{g}^{(d + 1)} = g^{(d + 1)} + θ (g^{(d + 1)} - g^{(d)})

ξ^{(d + 1)} = P_{A} (ξ^{(d)} - t α R \bar{g}^{(d + 1)})

Calculate the deviation value RD

R D = \frac{{‖ g^{(i + 1)} - g^{(i)} ‖}_{2}}{{‖ g^{(i + 1)} ‖}_{2}}

Until

R D \leq 10^{- 4}

, or

k

reaches the maximum number of iterations

Output recovery image g⁽ⁱ⁺¹⁾

where s is the iteration step of the original variable, t is the iteration step of the dyadic variable, and k is the number of iterations.

Then, there is $P_{A} (q) = \frac{q}{\max (1, {‖ q ‖}_{2}^{2})}$ .

2.3.3

Adaptive weighted total variation algorithm

The conventional TV term is based on the assumption that the reconstructed image is piecewise constantly distributed, and this approach can cause excessive smoothing of the reconstructed image edges. To alleviate the problem of over-smoothing of traditional TV edges, many researchers have studied weighted adaptive total variation [35]; therefore, in this study, an adaptive weighted total variation (AwTV) minimization image model was used for comparison with the proposed algorithm. $\min_{μ \geq 0} {‖ μ ‖}_{AwTV} subject to | p - A μ | \leq ε$ (30) where AwTV denotes the adaptive weight of the total variance of the reconstructed image. The ${‖ μ ‖}_{AwTV}$ is defined as follows: $\begin{array}{l} {‖ μ ‖}_{AwTV} \\ = \sum_{x, y, z} \sqrt{w_{x} {(μ_{x, y, z} - μ_{x - 1, y, z})}^{2} + w_{y} {(μ_{x, y, x} - μ_{x, y - 1, z})}^{2} + w_{z} {(μ_{x, y, z} - μ_{x, y, z - 1})}^{2}} \end{array}$ (31) $\begin{array}{l} w_{x} = \exp [- {(\frac{μ_{x, y, z} - μ_{x - 1, y, z}}{τ})}^{2}] \\ w_{y} = \exp [- {(\frac{μ_{x, y, z} - μ_{x, y - 1, z}}{τ})}^{2}] \\ w_{z} = \exp [- {(\frac{μ_{x, y, z} - μ_{x, y, z - 1}}{τ})}^{2}] \end{array}$ (32) where $w_{i}$ denotes the weight factor, and $τ$ is a scaling factor in the weights used to control the diffusion intensity of each iteration.

2.3.4

Fast Gradient Projection Algorithm

Beck and Teboulle proposed a fast gradient projection algorithm [34], which is derived as follows. The model expression to solve the TV-based denoising problem is given by Eq. (33). $\min_{x \in C} x - b_{F}^{2} + 2 λ TV (x)$ (33) TV is a nonsmooth regularization function that can be an isotropic ${TV}_{I}$ or anisotropic ${TV}_{l_{1}}$ .and $C$ belongs to a closed convex subset of $E \equiv R^{m \times n}$ . The nonsmoothness of the TV function leads to significant difficulties in solving Eq. (33). To address this challenge, Chambolle proposed a gradient-based algorithm for solving pairwise problems using a pairwise approach under unconstrained conditions. Thus, a constrained problem pair was constructed according to the proposed method. The following symbols are required.

• $P$ denotes the set of matrix-pairs $(p^{1}, p^{2})$ , which meet: $\begin{array}{l} {(p_{i, j}^{1})}^{2} + {(p_{i, j}^{2})}^{2} \leq 1, (0 \leq i \leq m, 0 \leq j \leq n) \\ | p_{i, n + 1}^{1} | \leq 1, 0 \leq i \leq m \\ | p_{m + 1, j}^{2} | \leq 1, 0 \leq j \leq n \end{array}$

• $L : ℝ^{(m - 1) \times n} \times ℝ^{m \times (n - 1)} \to ℝ^{m \times n}$ is linear operation, which is defined according to the following equation. $\begin{array}{l} L {(p^{1}, p^{2})}_{i, j} = (p_{i, j}^{1} - p_{i - 1 j}^{1}) + (p_{i, j}^{2} - p_{i, j - 1}^{2}) \\ 1 \leq i \leq m, 1 \leq j \leq n \end{array}$

• The operator $L^{T} : ℝ^{m \times n} \to ℝ^{(m - 1) \times n} \times ℝ^{m \times (n - 1)}$ adjacent to $L$ is calculated according to the following formula. $L^{T} (x) = (p^{1}, p^{2})$ where $p \in ℝ^{m \times n}$ is the matrix, whose defining equations are $\begin{array}{l} p_{i, j}^{1} = x_{i, j} - x_{i + 1, j}, i = 1, \dots, m - 1, j = 1, \dots, n \\ p_{i, j}^{2} = x_{i, j} - x_{i, j + 1}, i = 1, \dots, m, j = 1, \dots, n - 1 \end{array}$

P_C belongs to an orthogonal projection operator in set $C$ .

Therefore, when $C = B_{l, u}$ , the $P_{B_{l, u}}$ can be defined as: $P_{B_{l, u}} {(X)}_{i, j} = ｛ \begin{array}{l} \begin{array}{l} l x_{i j} < l \\ x_{i j} l \leq x_{i j} \leq u \\ u x_{i j} > u \end{array} \end{array}$

Based on the notation above, we obtain the dual problem in Eq. (33), and show the relationship between the primal and dual optimal solutions. To analyze this relationship, we propose the following hypothesis: We assume that $((p^{1}, p^{2})) \in P$ is the optimal solution to this problem. $\min_{(p^{1}, p^{2}) \in P} ｛ \begin{array}{l} h (p^{1}, p^{2}) \equiv - {‖ H_{C} (b - λ L (p^{1}, p^{2}) ‖}_{F}^{2} + \\ {| | b - λ L (p^{1}, p^{2}) | |}_{F}^{2} \end{array} ｝$ (34)

The $H_{C} (x)$ defined as: $H_{C} (x) = x - P_{C} (x), x \in ℝ^{m \times n}$ (35)

Because objective function $h$ in Eqs. (34) is continuously differentiable and its gradient is defined as $\nabla h (p^{1}, p^{2}) = - 2 λ L^{T} P_{C} (b - λ L (p^{1}, p^{2})),$ (36) $T (x, p^{1}, p^{2}) = T r (L {(p^{1}, p^{2})}^{T} x) .$ (37)

Therefore, problem Eq. (34) can be converted into the following equation. $\min_{x \in C} \max_{(p^{1}, p^{2}) \in P} ｛ X - b_{F}^{2} + 2 λ T r (L {(p^{1}, p^{2})}^{T} x) ｝$ (38)

Because the objective function is concave in $p^{1}, p^{2}$ and convex in $x$ , the maximum and minimum orders can be exchanged. Thus, the following equation is obtained: $max_{(p^{1}, p^{2}) \in P} \min_{X \in C} ｛ x - b_{F}^{2} + 2 λ T r (L {(p^{1}, p^{2})}^{T} x) ｝$ (39)

This can be rewritten as follows: $max_{(p^{1}, p^{2}) \in P} \min_{X \in C} ｛ X - {(b - λ L (p^{1}, p^{2}))}_{F}^{2} - b - λ L {(p^{1}, p^{2})}_{F}^{2} + b_{F}^{2} ｝,$ (40) where the optimal solution of the minimum problem is $x = P_{C} (b - λ L (p^{1}, p^{2})) .$ (41)

Subsequently, by substituting Eq. (41) into Eq. (40) and omitting the constant term, we obtain the following dual problem: $max_{(p^{1}, p^{2}) \in P} ｛ \begin{matrix} P_{C} (b - λ L (p^{1}, p^{2})) - (b - λ L {(p^{1}, p^{2})}_{F}^{2} \\ - b - λ L {(p^{1}, p^{2})}_{F}^{2} \end{matrix} ｝ .$ (42)

Our goal is to solve the dual problem in Eq. (34), whose gradient is expressed by Eq. (36). Therefore, the gradient projection algorithm, which is an effective method for solving the denoising problem, is presented. Because the norm of $x ϵ ℝ^{m \times n}$ is Frobenius norm. For $(p^{1}, p^{2}) ϵ ℝ^{(m - 1) \times n} \times ℝ^{m \times (n - 1)}$ , the norm is expressed as follows: $(p^{1}, p^{2}) = \sqrt{p^{1}_{F}^{2} + p^{2}_{F}^{2}},$ (43) $L h \leq λ^{2} .$ (44)

The projection onto set $P$ can be computed simply: For $(p^{1}, p^{2})$ , projection $P_{P} (p^{1}, p^{2})$ is given by $P_{P} (p^{1}, p^{2}) = (r^{1}, r^{2})$ . The equations of $r ϵ ℝ^{m \times n}$ are shown below: $\begin{array}{l} r_{i, j}^{1} = ｛ \begin{array}{l} \frac{p_{i, j}^{1}}{max ｛ 1, \sqrt{{(p_{i, j}^{1})}^{2} + {(p_{i, j}^{2})}^{2}} ｝} \\ \frac{p_{i, j}^{1}}{max ｛ 1, | p_{i, j}^{1} | + | p_{i, j}^{2} | ｝} \end{array} \\ r_{i, j}^{2} = ｛ \begin{array}{l} \frac{p_{i, j}^{2}}{max ｛ 1, \sqrt{{(p_{i, j}^{1})}^{2} + {(p_{i, j}^{2})}^{2}} ｝} \\ \frac{p_{i, j}^{2}}{max ｛ 1, | p_{i, j}^{1} | + | p_{i, j}^{2} | ｝} \end{array} \end{array}$ (45)

Substituting the objective function, gradient equation, and maximum Lipschitz constant into Eqs. (34) into the gradient projection algorithm introduced above, we obtain the following denoising algorithm:

2.3.5

OS-SART algorithm

In 1984, Anderson and Kak proposed an improved SART. The SART improves on the ART algorithms by calculating the error of all rays passing through a pixel at the same projection angle. The ART uses only one ray per iteration, whereas the SART corrects all the rays at the same angle. This is similar to lowering the noise introduced by the ART algorithm. Its formula is as follows: $x_{j}^{(k + 1)} = x_{j}^{k} + \frac{\sum_{p_{i \in} p_{φ}} ｛ λ_{k} \frac{p_{i} - \sum_{m = 0}^{M - 1} w_{i m} x_{i m}^{(k)}}{\sum_{m = 0}^{M - 1} w_{i m}} w_{i j} ｝}{\sum_{p_{i} \in p_{φ}} w_{i j}} .$ (46)

The OS-SART is a combination of ordered subsets and the SART algorithm, and its formula is as follows: $x_{j}^{(k + 1)} = x_{j}^{k} + \frac{\sum_{i \in S_{t}} ｛ λ_{k} \frac{p_{i} - \sum_{m = 0}^{M - 1} w_{i m} x_{i m}^{(k)}}{\sum_{m = 0}^{M - 1} w_{i m}} w_{i j} ｝}{\sum_{i \in S_{t}} w_{i j}},$ (47) where i denotes the i-th ray at a certain projection angle, j denotes the pixel value in the i-th ray, k denotes the number of iterations, and $λ$ is the relaxation factor.

2.3.6

Proposed algorithm

Based on the OS-SART iterative algorithm and first-order primal-dual algorithm with total variation denoising, we propose a new algorithm for sparse-view NCT 3D image reconstruction (OS-SART-PDTV). The iterative process of this algorithm contains two loops: the outer and inner loops are OS-SART for image reconstruction and PDTV for image denoising (Table 2).

OS-SART-PDTV algorithm

1. While the stop criterion is not met

2. Step1: OS-SART reconstruct the image

3. Initialization:

g_{0}

λ_{OS-SART}

N_{OS-SART}

λ_{PDTV}

and N_PDTV

4. For

n = 1

N_{OS-SART}

do;(n denotes the number of iterations)

g_{j}^{(k + 1)} = g_{j}^{k} + \frac{\sum_{i \in S_{t}} ｛ λ_{k} \frac{p_{i} - \sum_{m = 0}^{M - 1} w_{i m} x_{i m}^{(k)}}{\sum_{m = 0}^{M - 1} w_{i m}} w_{i j} ｝}{\sum_{i \in S_{t}} w_{i j}}

6. Non-negativity constraint;

g_{OS-SART}^{(n)} < 0, g_{OS-SART}^{(n)} = 0

7. Step2: PDTV

8. Set

s > 0, t > 0, θ \in [0, 1], input, u, α, R

9. Initialize

ξ^{(0)}, g^{(0)} = g_{OS-SART}

10. While the stop criterion is not met

11.

g^{(d + 1)} = {(H^{T} H + \frac{1}{s} I)}^{- 1} (H^{T} u + \frac{1}{s} g^{(d)} - α R^{T} ξ^{(d)})

12.

{\bar{g}}^{(d + 1)} = g^{(d + 1)} + θ (g^{(d + 1)} - g^{(d)})

13.

ξ^{(d + 1)} = P_{A} (ξ^{(d)} - t α R \bar{g}^{(d + 1)})

14. Calculate the deviation value RD,

R D = \frac{{‖ g^{(d + 1)} - g^{(d)} ‖}_{2}}{{‖ g^{(d + 1)} ‖}_{2}}

15. Until

R D \leq 10^{- 4}

, or

d

reaches the maximum number of iterations

16. Output recovery image

g^{(d + 1)}

17. end if;

18. Until the stopping criteria are reached.

19. Get the final reconstructed image g_OS-SART-PDTV

Experiment

3.1

Quantitative Evaluation Index

In this study, we used four image evaluation metrics to quantitatively analyze the mass of the reconstructed images for each algorithm. The accuracy of the image reconstruction was quantitatively analyzed using the correlation coefficient (CC) metric, which is defined as follows: $C C = \frac{\sum_{k = 1}^{M} (ϖ (y) - \bar{ϖ}) (ϖ_{t r u e} (y) - {\bar{ϖ}}_{t r u e})}{{｛ \sum_{k = 1}^{M} {(ϖ (y) - \bar{ϖ})}^{2} {\sum_{y = 1}^{M} (ϖ_{t r u e} (y) - {\bar{ϖ}}_{t r u e})}^{2} ｝}^{1 / 2}},$ (48) where ${\bar{ϖ}}_{t r u e}$ is the image to be reconstructed, M denotes the voxel number, and $ϖ (y)$ represents the reconstructed value at voxel $y$ . When the reconstructed image was the same as the image to be reconstructed, the CC value was 1.

The ((MSSIM) provides a comprehensive evaluation of the reconstructed image by comparing the differences between the image to be reconstructed and the reconstructed image in terms of brightness, contrast, and structure.

Brightness information: $l (ε, η) = \frac{2 u_{ε} u_{η} + c_{1}}{u_{ε}^{2} + u_{η}^{2} + c_{1}}$ .

Contrast information: $c (ε, η) = \frac{2 σ_{ε} σ_{η} + c_{2}}{σ_{ε}^{2} + σ_{η}^{2} + c_{2}}$ .

Structure Information: $s (ε, η) = \frac{ς_{ε η} + c_{3}}{ς_{ε} + ς_{η} + c_{2}}$ .

where $ε, η$ denote the standard image and the image to be evaluated, respectively; $u_{ε}$ and $u_{η}$ are the mean values of the image; $σ_{ε}$ and $σ_{η}$ are the standard deviations of the image; $ς_{ε η}$ is the covariance of the image; $C_{1}$ , $C_{2}$ and $C_{3}$ are constants. The mean structural similarity of an image is defined as $MSSIM = \frac{1}{M} \sum_{t = 1}^{M} {[l (ε_{t}, η_{t})]}^{a} \times {[c (ε_{t}, η_{t})]}^{β} \times {[s (ε_{t}, η_{t})]}^{γ} .$ (49)

RMSE is defined by the following equations. $RMSE = \sqrt{\frac{\sum_{n = 1}^{N} {({\hat{p}}_{n} - p_{n})}^{2}}{N}},$ (50)

${\hat{p}}_{n}$ represents a voxel in the image to be reconstructed, $p_{n}$ represents a voxel in the reconstructed sample, and $N$ indicates the number of voxels. A higher RMSE value indicates a larger error between the reconstructed images.

Universal quality image (UQI) is a widely used image evaluation index, which is defined as follows: $UQI = \frac{2 cov (\hat{μ}, μ)}{{\hat{ρ}}^{2} + ρ^{2}} \cdot \frac{2 \hat{μ} μ}{{\hat{μ}}^{2} + μ^{2}},$ (51) where cov denotes the covariance function, $\hat{μ}, μ$ are the means, and ${\hat{ρ}}^{2}, ρ^{2}$ are the variances of the reconstructed image and the image to be reconstructed, respectively. The UQI value ranges from 0 to 1 and increases with similarity.

3.2

Digital 3D Shepp–Logan Model Experiment

In the digital simulation experiment, we compared and analyzed the performance of five algorithms: FBP, OS-SART-TV, OS-SART-AwTV, OS-SART-FGPTV, and OS-SART-PDTV by reconstructing 3D Shepp–Logan model projection data. The dimensions of the 3D Shepp–Logan model were 256×256×256. An iterative algorithm must set reasonable iterative parameters to obtain an ideally reconstructed image. The following parameters were used to compare the three iterative algorithms: OS-SART-TV =2.7, Niter_OS-SART=80, $λ_{TV} = 200$ , and ${Niter}_{TV} = 300$ . For the OS-SART-AwTV algorithm =2.5, ${Niter}_{OS-SART} = 80$ , $λ_{AwTV} = 100$ , and ${Niter}_{TV} = 250$ . For the OS-SART-FGPTV algorithm, = 2.2, ${Niter}_{OS-SART} = 80$ , $λ_{PDTV} = 0.006$ , ${Niter}_{FGPTV} = 150$ . For the OS-SART-PDTV algorithm, = 2.2, ${Niter}_{OS-SART} = 80$ , $λ_{PDTV} = 0.0032$ , ${Niter}_{PDTV} = 300$ . To verify the denoising capability of the algorithm, the following parameters were set according to the noise model of the NCT system: These parameters include the photon incident flux 1×10⁵, m_ic=0, and $δ_{ic}^{2} = 10$ .

Fig.2 shows the relationship between the RMSE of the reconstructed images using the OS-SART-PDTV algorithm and the number of iterations for different sparse views (28 and 38 projection views), showing that the OS-SART-PDTV algorithm can converge to a stable solution after a certain number of iterations.

Fig. 2

RMSE versus iteration steps for the OS-SART-PDTV algorithm from different sparse views: (a) 28; (b) 38.

Fig. 3 shows the reconstruction of images with different numbers of projected views using the five algorithms. As illustrated in Fig. 3, the smoothness and sharpness of the images reconstructed using the five algorithms improved rapidly as the number of projection views increased. Compared to the other three iterative algorithms, the reconstructed images of the FBP algorithm showed severe bar artifacts. The reconstructed images from the four iterative algorithms are smoother and have sharper boundaries than those from the FBP algorithm. As shown in the enlarged ROI in Fig. 3, the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms can effectively reduce artifacts; however, part of the edge structure is not exactly reconstructed. In contrast, the novel algorithm can not only maintain good reconstruction performance but also reconstruct finer structures.

Fig. 3

Sliced images at z=130 for four algorithms reconstructing reconstructed images with different numbers of projection views.

In addition, we compared and analyzed the error images of the reconstructed images using different algorithms that reflected the difference between the pixel values of the reconstructed and reference images. According to the error image in Fig. 4, the image reconstructed using the FBP algorithm lost many features and produced many artifacts, indicating that the algorithm is not suitable for sparse-view NCT 3D reconstruction. The other four iterative algorithms performed better at suppressing artifacts and reconstructing detailed structural information. The OS-SART-PDTV algorithm has much less detail loss in the error image than the other algorithms.

Fig. 4

Slice of error images at z=130

We plotted the profile of the reconstructed images for each algorithm in Fig. 5 to further evaluate the performance of each algorithm. The number of projection views from left to right was 28, 38, and 38 +, respectively. Here, 38+ indicates that the noise model was added to the 3D Shepp-Logan model. As shown in Fig. 5, there was a large fluctuation in the FBP profile line with the largest deviation from the true pixel value. Although the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms performed better than the FBP algorithm, there was still a gap between their amplitudes and true pixel values. The profile of OS-SART-PDTV matches the reference value best and is closest to the true value. This also indicates that the OS-SART-PDTV algorithm is well-suited for sparse-view NCT 3D reconstruction.

Fig. 5

(Color online) Horizontal profiles of the 3D Shepp–Logan model at z =130^th slices. (a) 28; (b) 38; (c) 38+.

As shown in Fig. 6. We used four image evaluation metrics to quantitatively analyze the reconstructed images from each algorithm. First, the evaluation indices of the images reconstructed using noiseless projection data were analyzed. From the CC, MSSIM, and UQI histograms shown in Figure 6, the values of the three metrics for the reconstructed images from each reconstruction algorithm gradually increased as the number of projection views increased. These three metrics exhibit similar characteristics.

Fig. 6

(Color online) Histogram of reconstructed image evaluation metrics for the 3D Shepp-Logan model. (a) CC; (b) MSSIM; (c) RMSE; (d) UQI.

According to the RMSE histogram in Fig. 6c, the reconstructed images of OS-SART-PDTV had the smallest RMSE values compared to the other algorithms for the same projection views. The RMSE of the reconstructed image for each algorithm gradually decreased as the number of projection views increased. The RMSE value of the images reconstructed using the OS-SART-PDTV algorithm at 28 projection views was 0.03352, indicating that the reconstructed image was closest to the real image and that the quality of the reconstructed images was higher. As shown in the UQI histogram in Fig. 6d, the UQI values of the reconstructed images for each algorithm gradually increased as the number of projected views increased. The trend of the UQI was the opposite to that of the RMSE. The reconstructed images from the OS-SART-PDTV algorithm had the highest UQI values for the same number of projection views, indicating that it outperformed the other algorithms. The UQI value of the image reconstructed using the OS-SART-PDTV algorithm was 0.99039 for the 28 projection views. We then quantitatively analyzed four evaluation metrics for the images reconstructed using data containing noisy projections. From the CC, MSSIM, and UQI histograms, it can be seen that the OS-SART-PDTV algorithm has the largest evaluation metrics for these three reconstructed images when the same number of noise-containing projections are reconstructed. According to the RMSE histogram, the reconstructed image from the OS-SART-PDTV algorithm had the lowest RMSE value. Therefore, based on the quantitative and visual analyses of each algorithm, the OS-SART-PDTV algorithm has the highest-quality reconstructed images with the same number of projection views.

3.3

Digital Head Model Experiment

We further analyzed the performance of these five algorithms by reconstructing a digital head model. The reconstructed images from the five algorithms used for reconstructing different projection view numbers are shown in Fig.7. Based on the reconstructed images, the smoothness and sharpness of the reconstructed images of all five algorithms improve significantly as the number of projected views increases. Compared to the other four algorithms, the reconstructed images of the FBP algorithm show severe bar artifacts. The reconstructed images from these four iterative algorithms were smoother and had cleaner image boundaries than those from the FBP algorithm. As shown in the enlarged ROI in Fig.7, the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms are effective in reducing artifacts, but the edge structure of the image is not accurately reconstructed. In contrast, the proposed OS-SART-PDTV algorithm can not only reduce artifacts, but also reconstruct more detailed structural information.

Fig. 7

(Color online) Sliced images at z=70 for four algorithms reconstructing reconstructed images with different number of projection views.

To further evaluate the performance of the algorithm, the contours of the reconstructed images for each algorithm are plotted in Fig. 8. As shown in Fig. 8, the numbers of projected views from left to right are 28, 38, and 38+, where 38+ indicates that noise has been added to the projection data. As shown in Fig.8, there was a large fluctuation in the FBP profile, which had the largest gap from the true pixel values. Although the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms performed better than the FBP algorithm, there was still a large deviation between the actual and true pixel values. The profile of OS-SART-PDTV matches the reference value best and is closest to the true value. This also indicates that the OS-SART-PDTV algorithm is well-suited for sparse-view NCT 3D reconstruction.

Fig. 8

(Color online) Horizontal profiles of the digital head model at z=250^th slices. (a) 28; (b) 38; (c) 38+.

According to the RMSE histogram of the reconstructed images from each algorithm, the OS-SART-PDTV algorithm had the smallest RMSE of the reconstructed images. As the number of projected views gradually increased, the RMSE of the reconstructed image for each algorithm gradually decreased. When the number of projected views is 28, the RMSE value of the OS-SART-PDTV algorithm-reconstructed image is 0.02404, which indicates that the reconstructed image of this algorithm is the closest to the real image, and the quality of the reconstructed image is higher.

The UQI histogram in Fig. 9d shows that the UQI values of the reconstructed images of each algorithm increased as the number of projected views increased. The UQI trend was opposite to that of the RMSE. The OS-SART-PDTV algorithm reconstructs the image with the largest UQI value for the same number of projected views, indicating that this algorithm outperforms the other three iterative algorithms. When the number of projected views was 28, the UQI of the image reconstructed by PDTV algorithm was 0.98594. We then quantitatively evaluated images containing noisy projections reconstructed using the four algorithms. The OS-SART-PDTV algorithm reconstructs the image with the largest CC, MSSIM, and UQI values for the same number of projected views. According to the RMSE histogram, the RMSE value of the reconstructed image obtained using OS-SART-PDTV was the smallest. Therefore, by quantitative and visual analyses of each algorithm, OS-SART-PDTV exhibited superior performance in sparse view reconstruction.

Fig. 9

(Color online) Histogram of reconstructed image evaluation metrics for the digital head model. (a) CC; (b) MSSIM; (c)RMSE; (d)UQI.

To analyze the relationship between the OS-SART-PDTV algorithm and regularization parameters, we present the reconstructed images with different regularization parameters in Fig.10. Based on the reconstructed images, an optimal image was obtained by adjusting the regularization parameters when reconstructing the head model. According to the reconstructed images with different regularization parameters, it can be seen that the reconstructed image quality reaches the best when the regularization parameters reach a certain value as the reconstruction parameters gradually increase. As the regularization parameter continued to increase, the quality of the reconstructed image decreased.

Fig. 10

(Color online) Reconstructed images with different regularization parameters for the OS-SART-PDTV algorithm 38+

3.4

Real neutron projection experiment

To further validate the performance of this novel algorithm in practical applications, the performance of the algorithm was verified using projection data from an NCT-based clock model (Fig.11) provided by Burkhard Schillinger of the Technical University of Munich, Germany [35]. Schillinger et al. provided 201 neutron projection images captured at equal angles in the range of 0 °to 180°. They also provided two dark field and two open field background images. The background images were obtained using the CCD camera with the neutron beam turned off, which contained readout noise as well as electron noise caused by dark currents. Open-field images were obtained using a CDD camera with neutrons turned on and without the samples.

Fig. 11

(Color online) Clock model

As shown in Fig.12, the artifacts of the reconstructed images by the FBP algorithm gradually increase as the number of projected views decreases. When the number of projected views is 45, the four iterative algorithms outperform the FBP algorithm in suppressing artifacts and reconstructing fine structure information. Although the OS-SART-TV, reconstructs images with fewer artifacts, the edges of the reconstructed images are too smooth, resulting in the loss of detailed structural information of the reconstruction images. Despite OS-SART-AwTV and OS-SART-FGPTV algorithms can also reconstruct a lot of image details, their reconstruction results are poor compared with the OS-SART-PDTV algorithm that reconstructs the same number of sparse views. Visually, the OS-SART-PDTV algorithm outperforms other algorithms in terms of noise removal, artifact suppression, and retention of detailed structural information.

Fig. 12

(Color online) Sliced images at z=201 for four algorithms reconstructing reconstructed images with different number of projection views.

The performances of the three iterative algorithms were quantitatively analyzed using four image-evaluation metrics, as shown in Fig.13. In Fig.13, the RMSE of the image reconstructed by the OS-SART-PDTV algorithm is the smallest, indicating the minimal error between the reference image and the reconstructed image, whereas the image reconstructed by the OS-SART-PDTV algorithm has the highest CC, MSSIM, and UQI values, indicating that it is most similar to the reference image. In summary, the OS-SART-PDTV algorithm performed well in denoising, suppressing artifacts, and reconstructing fine structures when reconstructing the projection data of sparse NCT scans, indicating that this novel algorithm can be applied to the 3D reconstruction of sparse-view NCT.

Fig. 13

(Color online) Histogram of reconstructed image evaluation metrics for the clock model. (a) CC; (b) MSSIM; (c)RMSE; (d)UQI.

Figure 14 presents the profiles of the reconstructed images. According to the profile, the OS-SART-PDTV algorithm has an advantage over the other algorithms, and its profile is the closest to the reference.

Fig. 14

(Color online) The profile of reconstructed images by FBP, OS-SART-TV, OS-SART-AwTV, OS-SART-FGPTV, and OS-SART-PDTV at z=201. The number of projection views is (a) 45; (b) 90; (c) 135.

Discussion

Herein, we proposed a 3D reconstruction algorithm for sparse-view NCT. Because the analytical algorithm (FBP) has inherent defects that make it unsuitable for sparse-view reconstruction, we used an iterative algorithm applied to the 3D reconstruction of sparse-view NCT. Furthermore, the iterative algorithm can embed a priori information to regularize the image. Therefore, the TV-based iterative algorithm is widely used in sparse-view CT image reconstruction because it is superior to other algorithms in artifact suppression and denoising. By reconstructing the projection data of the 3D Shepp-Logan, head, and clock models, the OS-SART-PDTV algorithm outperformed the other algorithms in both quantitative and visual analyses. Therefore, OS-SART-PDTV is suitable for sparse-view NCT. In addition, deep learning, an important research direction in the field of image reconstruction, has led to many significant research results in sparse-view reconstruction, such as in [36], using short C-arm CT. scans. This study proposed the DRONE model, which uses a codec network in the embedding module to extract depth features in the data and image domains, combined with the Wasserstein distance generation adversarial network to maintain details and features in the image domain. It combines data residual and image residual networks in the refinement module to restore fine structure features from the output of the embedding module. According to the CS iterative reconstruction model, the depth prior of the data and image domains is regularized, thus ensuring the robustness of the DRONE network and making the DRONE model good at retaining image edge information, recovering image features, and reconstructing accurate image information in practical applications.

Generally, in practical NCT scanning systems, the projected image is affected by photon statistical noise and electrical noise. Therefore, several iterative reconstruction algorithms cannot achieve satisfactory results in NCT scanning systems. To further test the performance of our algorithm, we reconstructed the projection data obtained from the NCT scan-clock model using the OS-SART-PDTV algorithm. Based on the results of the reconstructed neutron projection data, OS-SART-PDTV exhibited a good performance in denoising, suppressing artifacts, and reconstructing fine structures. According to the results of the reconstruction of the 3D Shepp-Logan model, the OS-SART-PDTV algorithm converged monotonically to a stable solution (Fig. 3).

Iterative reconstruction algorithms generally need to optimize several parameters to obtain a stable solution. This is a common problem encountered by all iterative reconstruction algorithms. The OS-SART-PDTV algorithm must optimize the following parameters: $λ_{OS-SART}$ , $N_{OS-SART}$ , $λ_{PDTV}$ and $N_{PDTV}$ . We found that the quality of the reconstructed images was significantly improved when $λ_{OS-SART}$ varied within a small range. In other words, the OS-SART–PDTV algorithm was more sensitive to this parameter. $N_{OS-SART}$ is usually set in the range of 50–100 because the improvement rate of the image quality gradually decreases as the number of iterations increases, and after reaching a certain number of iterations, the image quality no longer improves. According to our experimental results, the OS-SART-PDTV algorithm obtained a better reconstructed image after 50 iterations. In addition, the two parameters $λ_{PDTV}$ and $N_{PDTV}$ must be tuned several times to obtain high-quality reconstructed images. When reconstructing the projection data of a new sample, the two parameters $λ_{PDTV}$ and $N_{PDTV}$ must be readjusted.

Conclusion

NCT has proven to be a highly effective noninvasive measurement method and has been widely applied in nuclear engineering, thermodynamics, and other fields. However, because of the complexity of the NCT application environment, obtaining complete projection data is challenging. The selection of an efficient sparse-view reconstruction algorithm is directly related to the system imaging quality and reliability of the measurement results. To this end, we propose the OS-SART-PDTV algorithm for sparse-view NCT 3D. The OS-SART-PDTV algorithm is divided into two steps: OS-SART for image reconstruction, and PDTV for solving the total variance using a first-order primal-dual algorithm. Based on the reconstruction results of the projection data obtained from the NCT scan clock model, the OS-SART-PDTV algorithm has a significant advantage over the other algorithms in preserving edge structure, denoising, and suppressing artifacts. Therefore, the OS-SART-PDTV algorithm has great potential for application in NCT 3D reconstruction.

Reference

F. Tian, C R. Geng, X B. Tang et al.,

Analysis of influencing factors on the method for determining boron concentration and dose through dual prompt gamma detection

. Nucl. Sci. Tech. 32(4), 35 (2021). doi: 10.1007/s41365-021-00873-3