Sparse-View Neutron CT Reconstruction Using a Modified Weighted Total Difference Minimization Method

Abstract

Indirect neutron imaging is an effective method for nondestructive testing of spent nuclear fuel elements. Considering the difficulty of obtaining experimental data in a high-radiation environment and the characteristic of high noise of neutron images, it is difficult to use the traditional FBP algorithm to recover the complete information of the sample based on the limited projection data. Therefore, it is necessary to develop the sparse-view CT reconstruction algorithm for indirect neutron imaging. In order to improve the quality of the reconstruction image, an iterative reconstruction method combining SIRT, MRP, and WTDM regularization is proposed. The reconstruction results obtained by using the proposed method on simulated data and actual neutron projection data are compared with the results of four other algorithms (FBP, SIRT, SIRT-TV, and SIRT-WTDM). The experimental results show that the SIRT-MWTDM algorithm has great advantages in both objective evaluation index and subjective observation in the reconstruction image of simulated data and neutron projection data.

1. Introduction

With the rapid growth of the world economy, non-renewable resources such as oil, natural gas, and coal are constantly consumed, and the problem of the energy crisis has become more and more obvious [1]. To seek the way of sustainable development, nuclear power plants have become a hotspot of people’s attention in many countries. At present, nuclear power plants have become one of the most important components of modern energy [2,3].

Nuclear safety is the lifeline for the healthy development of nuclear power, which directly affects people’s normal life, production, and social activities. Especially after the Three Mile Island nuclear power plant leakage accident in the United States in 1979, the Chernobyl nuclear power plant accident in the former Soviet Union in 1986, and the Fukushima nuclear power plant accident in Japan in 2011, the public has paid great attention to the safety of nuclear power plants [4,5,6]. The nuclear fuel element is the core component of a reactor. The nuclear fuel pellet located in the cladding releases energy through nuclear fission. The interaction between the pellet and the cladding will affect the safety of the fuel element. If the pellet is damaged, the debris will enter the gap between the pellet and the cladding, which may cause close contact between the pellet and the cladding to form a hotspot. The hotspot will result in the local temperature of the cladding being too high and the rupturing of the cladding, which can lead to nuclear leakage. Therefore, the safety inspection is one of the important means to ensure the safe operation of nuclear power plants.

Nondestructive testing technology has the advantages of being nondestructive and rapid and reflecting the overall structure information of the tested samples [7]. In recent years, with the development of electronics and information science, X-ray and neutron imaging technologies have become more and more mature. Advanced detection imaging equipment and image processing software have been developed and applied continuously and can realize three-dimensional tomography and visualization of the tested samples [8,9]. However, the spent nuclear fuel element is a strong emitter of X-rays and gamma rays, and it is impossible to analyze the irradiated object by X-ray computed tomography (CT). Neutron CT is useful for nondestructive testing of spent nuclear fuel elements due to the suitable neutron attenuation properties of uranium [10].

When scanning in the range of 180 or 360°, if the detected object is sampled by a high-density projection, the traditional filtered back-projection (FBP) algorithm will be used to reconstruct high-quality images. When the number of projections does not satisfy the complete sampling condition required by Shannon’s sampling theorem for reconstruction, it is called sparse-view CT reconstruction [11]. The iterative reconstruction algorithm is mainly used for image reconstruction of sparse-view projections, such as iterative filtered back-projection algorithm, algebra reconstruction technique (ART) [12], simultaneous iterative reconstruction technique (SIRT) [13], and simultaneous algebraic reconstruction technique (SART) [14]. With the development of compressed sensing theory, algebraic iterative reconstruction theory based on total variation (TV) has great potential in the reconstruction of sparse-view projections [15,16]. Rudin proposed an ROF denoising model for the image total variation, which used the total variation as a regular constraint for image denoising so that the edge and detail information of the image could be well preserved while the noise is reduced [17]. Sidky et al. proposed TV-POCS and ASD-POCS algorithms, which have played a good role in the aspect of sparse-view reconstruction [18,19]. Several methods based on TV and its variants, including edge-preserving TV [20], anisotropic TV [21,22], weighted TV [23,24], and directional TV [25], were proposed to improve the performance of the algorithm.

The regularized iterative reconstruction algorithm based on total difference minimization (TDM) combined with soft-threshold filtering has achieved a good effect for sparse-view CT reconstruction [26]. However, the total difference (TD) regularization only considers the sparsity of the image gradient in the horizontal and vertical directions, and the protection of the structure information is a little insufficient at the edge of the object. Shu and Ahuja considered a new regularization measure method based on compressed sampling, which constrains both the sparsity and continuity of the gradient and is called weighted total difference minimization (WTDM) [27]. According to their research, the hybrid compression sampling algorithm can obtain a high-precision image, and the edges in the image can be well restored. The weighted total difference minimization not only considers the sparsity constraint of the image gradient, but also considers the continuity constraint of the image gradient, which can effectively protect the gradient information in different directions of the image edge. In addition, the neutron projected data contain isolated noise, which affects the quality of a reconstructed image. Median root prior (MRP) regularization can suppress local non-monotonic structures and eliminate the influence of isolated noise on image quality [28,29]. Inspired by MRP and WTDM regularization algorithms, this paper proposes a method that can simultaneously remove the isolated noise and effectively protect the gradient information in different directions. To make the soft threshold filter function adaptively adjust the threshold according to the image grayscale, the threshold rule is modified in this paper.

The process of the proposed algorithm (simultaneous iterative reconstruction technique—modified weighted total difference minimization (SIRT-MWTDM)) in this paper is that the image reconstructed by using the FBP algorithm is used as the initial image of the SIRT-MWTDM algorithm. After each SIRT iteration, MRP and WTDM regularization operations are performed immediately. The loop is repeated until the stop condition is met. The main contribution of this paper is the improvement of the WTDM soft threshold selection rules, and MRP and WTDM regularization are used in our algorithm. The SIRT-MWTDM algorithm is applied to the sparse angle reconstruction from the high-noise projections. The results of the simulation and neutron experiment show that our algorithm is superior to other algorithms in suppressing noise and improving image clarity. Moreover, it allows high-quality reconstructed images to be provided to the researchers of nuclear fuel elements in a reactor, which is very meaningful.

2. Methods

2.1. Weighted Total Difference Minimization

Because the projection data used in sparse-view CT reconstruction are incomplete, it is an ill-posed inverse problem. At present, a large number of documents realize the solution to such ill-posed problems by introducing the regularization method into the process of reconstruction. Finally, a small amount of projection data can be used to reconstruct high-quality CT images. The regularization method is used to solve this kind of ill-posed problem, which can be equivalent to the solution process of the optimization problem as shown in Equation (1).

$$f=\underset{f}{\arg \min }{\Vert Af-p\Vert }_2^2+\lambda \cdot R(f)$$

(1)

where $f$ is the image to be restored, $A$ is the projection matrix, $p$ is the projected image which is captured, $\lambda$ is a regularization parameter, ${\Vert Af-p\Vert }_2^2$ is a fidelity item, and $R(f)$ is the regular term of the iterative reconstruction algorithm.

According to the literature [30], the weighted total difference minimization considers the sparsity constraint of the image gradient and the continuity constraint of the image gradient, which can effectively protect the gradient information in different directions at the edge of the image. Therefore, we take WTDM as a regularization constraint for the iterative reconstruction of sparse-view neutron images. According to WTD theory, the WTD of an image $f$ can be defined according to Equation (2):

$$\mathrm{WTD}(f)={\Vert G_{\mathrm{x}}{f}_{i,j}\Vert }_1+{\Vert G_{\mathrm{y}}{f}_{i,j}\Vert }_1+\phi ({\Vert G_{xy}{f}_{i,j}\Vert }_1+{\Vert G_{\mathrm{yx}}{f}_{i,j}\Vert }_1)$$

(2)

where

$$\begin{array}{l}{\Vert G_{\mathrm{x}}{f}_{i,j}\Vert }_1={\Vert f_{i+1,j}-f_{i,j}\Vert }_1\\ {\Vert G_{\mathrm{y}}{f}_{i,j}\Vert }_1={\Vert f_{i,j+1}-f_{i,j}\Vert }_1\\ {\Vert G_{\mathrm{yx}}{f}_{i,j}\Vert }_1={\Vert f_{i,j+1}-f_{i+1,j}\Vert }_1\\ {\Vert G_{\mathrm{xy}}{f}_{i,j}\Vert }_1={\Vert f_{i+1,j+1}-f_{i,j}\Vert }_1\end{array}$$

(3)

where $f_{i,j}$ represents the pixel value in the coordinates $(i,j)$ of the reconstructed image. $\phi$ is the weight value between gradient continuity and gradient sparsity constraints. When $\phi =0$, it only constrains the sparsity of the image gradient. When $\phi =1$, the same degree of constraint for gradient sparsity and continuity can be achieved. If we substitute the regularization term $\mathrm{WTD}(f)$ into the iterative reconstruction Equation (1), the image iterative reconstruction model is transformed into Equation (4).

$$f=\arg \underset{f}{\min }{\Vert Af-g\Vert }_2^2+\lambda \cdot {\displaystyle \sum _i{\displaystyle \sum _j\left[\left|f_{i+1,j}-f_{i,j}\right|+\left|f_{i,j+1}-f_{i,j}\right|+\phi (\left|f_{i+1,j+1}-f_{i,j}\right|+\left|f_{i,j+1}-f_{i+1,j}\right|)\right]}}$$

(4)

Therefore, the solution of the weighted total difference minimization CT reconstruction model is equivalent to the solution of Equation (5). There are many methods for solving Equation (5), such as the soft threshold filtering method and the split Bregman method [31]. The convergence and effectiveness of the soft threshold filtering method have been proved in theory, and the soft threshold filtering method has been successfully applied to the field of CT reconstruction [26]. Therefore, the soft threshold filtering method is also used in this paper to solve Equation (5). In the framework of soft threshold filtering, a pseudoinverse equation is constructed to solve the reconstruction model based on weighted total difference minimization.

$$f_{i,j}^{n+1}=\frac{1}{4+4\cdot \phi }\times \{\begin{array}{l}F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i+1,j}^{n+1})+F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i,j+1}^{n+1})+F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i,j-1}^{n+1})+F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i-1,j}^{n+1})\\ +\phi \cdot \left[F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i+1,j+1}^{n+1})+F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i+1,j-1}^{n+1})+F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i-1,j-1}^{n+1})+F(\omega ,{\tilde{f}}_{i,j}^{n+1},{\tilde{f}}_{i-1,j+1}^{n+1})\right]\end{array}\}$$

(5)

where ${\tilde{f}}_{i,j}^{n+1}$ is the gray value at the coordinate position $(i,j)$ of the reconstructed image.

$$F(\omega ,y,z)=\{\begin{array}{l}(y+z)/2,\begin{array}{cc}& \end{array}\left|y-z\right|<\omega \\ y-\omega /2,\begin{array}{cc}& \end{array}(y-z)\ge \omega \\ y+\omega /2,\begin{array}{cc}& \end{array}(y-z)\le -\omega \end{array}$$

(6)

where y and z are input parameters of the soft threshold function F.

To make the $\omega$ value change adaptively according to the local image, we redefine $\omega$ to replace the original fixed value, as shown in Equation (7).

$$\omega =\mathrm{abs}(f-\mathrm{Med}(f))/k$$

(7)

where $k(k>0)$ can be set according to the actual effect of the reconstructed image, which is mainly used to adjust the value of the threshold $\omega$.

2.2. Median Root Prior

The median filter is used as the energy function in the MRP algorithm so that the MRP algorithm can suppress isolated noise. The neutron projection images often contain isolated noise, which affects the quality of reconstructed images. Hence, the MRP algorithm is introduced into the reconstruction process to realize the local monotone constraint of the image and to suppress the influence of isolated noise. The MRP algorithm formula is shown in Equation (8), and the noise suppression formula of the reconstructed image is shown in Equation (9):

$${\eta }_{\mathrm{MRP}}^{(n)}(i,j)=\frac{1}{1+\gamma \frac{f^{(n)}(i,j)-\mathrm{Med}(f^{(n)}(i,j))}{\mathrm{Med}(f^{(n)}(i,j))+\epsilon }}$$

(8)

$$f_{\mathrm{MRP}}^{(n)}(i,j)={\eta }_{\mathrm{MRP}}^{(n)}(i,j)\times f^{(n)}(i,j)$$

(9)

where $n$ is the iteration number of the reconstruction algorithm. $\mathrm{Med}(f^{(n)}(i,j))$ is the median filtering algorithm by $3\times 3$ or $5\times 5$ pixel neighborhood. $f^{(n)}(i,j)$ is the reconstructed image after $n$ iterations. $f_{\mathrm{MRP}}^{(n)}(i,j)$ is the image of suppressed noise. $\epsilon$ is a perturbation parameter that prevents the denominator from being infinite, which is 10⁻⁶ in this paper. $\gamma$ is the noise suppression intensity parameter of the MRP algorithm with the value ranging from 0 to 1. If a reconstructed image contains locally nonmonotonic structures such as isolated noise spikes, the MRP factor is not equal to 1 because of the inequality between the pixel value and its median value. Then, the factor will modify the pixel value, and these structures will be suppressed. According to our research, if $\gamma =0$, the MRP is not used to suppress noise. A large $\gamma$ corresponds to a strong weight to suppress isolated noise.

2.3. The SIRT-MWTDM Iterative Algorithm

The SIRT algorithm adopts the point-by-point update mode, and the contribution value of rays from all projection angles to the current point should be calculated for updating each point. The update of each image point is obtained by averaging the contribution value of all rays passing through the point. It has more advantages in suppressing statistical noise than ART and SART. The iterative equation of the SIRT algorithm is shown in Equation (10).

$$f_j{}^{(n+1)}=f_j{}^{(n)}+{\lambda }_{\mathrm{SIRT}}^{(n)}\cdot {\displaystyle {\sum }_{i=1}^N[a_{ij}(p_i-{\displaystyle {\sum }_{m=1}^M{a}_{im}{f}_m})/{\displaystyle {\sum }_{m=1}^M{a}_{im}}]}/{\displaystyle {\sum }_{i=1}^N{a}_{ij}}$$

(10)

where $n$ is the iteration number of the SIRT algorithm and $f_j$ is the pixel value of the reconstructed image. Each ray beam is regarded as an infinitely narrow straight line, and the length of the ray passing through the reconstructed pixel is defined as $a_{ij}$. Namely, $a_{ij}$ is an element value of the system matrix. $1\le i\le N$, $1\le j\le M$, and $N$ and $M$ are the width and height of the sinogram. ${\lambda }_{\mathrm{SIRT}}$ is the relaxation factor ( $0<{\lambda }_{\mathrm{SIRT}}<2$).

Before the SIRT-MWTDM algorithm is used for iterative reconstruction, we use the FBP algorithm to obtain a reconstruction image $f^{(0)}$ from sparse angle projections. $f^{(0)}$ is used as the initial image of the SIRT-MWTDM algorithm. Then, the SIRT-MWTDM algorithm starts the process of iterative reconstruction according to the parameters that have been set. Each complete iteration process of the SIRT-MWTDM algorithm includes one SIRT iteration, one MRP, and one WTDM. When the stop condition is met, we can obtain the reconstructed image. In the initial stage of the SIRT-MWTDM algorithm, the estimated solution is far from the optimal solution, and a large ${\lambda }_{\mathrm{SIRT}}$ value can speed up the optimization process. As the iteration number increases, the estimated solution becomes closer and closer to the optimal solution. As the number of iterations increases, the smaller ${\lambda }_{\mathrm{SIRT}}$, the higher quality of the convergence. Therefore, we introduce the parameter ${\lambda }_{red}$ of attenuation rate to our algorithm. To more clearly illustrate the SIRT-MWTDM algorithm, the organizational scheme of the SIRT-MWTDM algorithm is shown in Figure 1 and Algorithm 1.

Algorithm 1. The organizational scheme of the SIRT-MWTDM

Initialization: ${\lambda }_{\mathrm{SIRT}}$; ${\lambda }_{red}$; $\phi$; $\gamma$; $N_{\mathrm{iter}}$; $k$;

FBP reconstruction:

$f^{(0)}=\mathrm{FBP}(p)$

for $n=1$ to $N_{\mathrm{iter}}$ do

SIRT updating:

${\lambda }_{\mathrm{SIRT}}^{(n)}={\lambda }_{\mathrm{SIRT}}^{(n-1)}\times {\lambda }_{red}$

$f^{(n)}=\mathrm{SIRT}(f^{(n-1)})$

Non-negativity constraint:

If $f^{(n)}<0$, $f^{(n)}=0$

MRP:

Calculate ${\lambda }_{\mathrm{MRP}}^{(n)}$ about image $f^{(n)}$ according to Equation (8)

Calculate $f_{\mathrm{MRP}}^{(n)}$ according to Equation (9)

WTDM:

Calculate adaptive threshold $\omega$ according to Equation (7)

Solve equation (4) according to Equations (5) and (6)

Image updating and next loop

end

3. Experiment

3.1. Quantitative Evaluation Index

To quantitatively compare the reconstructed effects of the different reconstruction algorithms from sparse-view projections, three indexes, namely structural similarity (SSIM), peak signal to noise ratio (PSNR), and root mean squared error (RMSE), are used for quantitative comparative analysis [32]. The equations of the three indexes are shown in Equations (11)–(13).

$$SSIM(f_A,f_B)=\frac{(2{\mu }_{f_A}{\mu }_{f_B}+c_1)(2{\sigma }_{f_A{f}_B}+c_2)}{({\mu }_{f_A}^2+{\mu }_{f_B}^2+c_1)({\sigma }_{f_A}^2+{\sigma }_{f_B}^2+c_2)}$$

(11)

where $f_A$ is the reference image and $f_B$ is the reconstructed image. ${\mu }_{f_A}$ and ${\mu }_{f_B}$ are the mean values of $f_A$ and $f_B$. ${\sigma }_{f_A}$ and ${\sigma }_{f_B}$ are the standard deviations of $f_A$ and $f_B$. ${\sigma }_{f_A{f}_B}$ is the covariance between $f_A$ and $f_B$. $c_1$ and $c_2$ are constants with a small value, which are used to maintain stability. In this paper, we set $c_1=c_2=0.001$. SSIM is used to measure the similarity between two images. A larger SSIM value indicates a higher similarity and a smaller structural difference between the two images.

$$PSNR=10{\log }_{10}\frac{{(\max (f_{\mathrm{A}}))}^2}{\frac{1}{M\times N}{\displaystyle {\sum }_{i=1}^M{\displaystyle {\sum }_{j=1}^N{(f_{\mathrm{A}}-f_{\mathrm{B}})}^2}}}$$

(12)

where $M$ and $N$ are the width and height of image $f_A$ and $f_B$. PSNR is used to measure the quality of the reconstructed image. If the PSNR value is larger, the quality of the reconstructed image is higher.

$$RMSE=\sqrt{\frac{1}{M\times N}{\displaystyle {\sum }_{i=1}^M{\displaystyle {\sum }_{j=1}^N{(f_{\mathrm{A}}-f_{\mathrm{B}})}^2}}}$$

(13)

RMSE is used to evaluate the quality of the reconstructed image and to measure the deviation between the reconstructed image and the reference image. A small RMSE value indicates that the difference between the two images is small.

3.2. Simulation Experiment

Due to the influence of neutron beam quality and hardware inherent noise, the quality of the reconstructed image is affected by noise in the collected projections. Therefore, by adding Poisson noise into the sinogram of the simulated models, we can not only verify the effectiveness of the sparse-view reconstruction algorithm, but also test the robustness of the reconstruction algorithm to the noise. We set the number of photons according to Beer’s law to add noise to the simulation images. The number of quanta is set $5\times {10}^4$ for Poisson noise. In the simulation experiment, the projected images are uniformly collected within the range $[0^{\mathrm{o}},{180}^{\mathrm{o}}]$ with equal intervals. The number of projections collected is 30 and 60, and the corresponding sampling angle intervals are 6^o and 3^o. We collected three sets of simulation data, namely 30 and 60 noiseless projections and 60 noisy projections. Because the neutron beams are parallel, there is no need to set the distance between the neutron source and the detector. Both models are $512\times 512$ in size as shown in Figure 2. The number of pixels of the detector is 512, and the pixel unit is mm. Figure 2a is the Shepp–Logan phantom. Figure 2b is a simulated circular phantom based on the structure of the nuclear fuel element. We provide the attenuation coefficients of each part in the Shepp–Logan and circular phantoms.

3.2.1. Shepp–Logan Phantom

Four reconstruction algorithms, FBP, SIRT, SIRT-TV, and SIRT-WTDM, are used for comparison with the reconstructed results obtained by using the SIRT-MWTDM algorithm in this paper. The calculation process of the SIRT-TV algorithm includes one SIRT and one TV. We can obtain the reconstructed image when the set number of iterations is reached. To make a reasonable and fair comparison of the four algorithms, we set the parameters of the algorithms as follows: ${\mathit{\lambda }}_{\mathrm{SIRT}}=1.8$, ${\mathit{\lambda }}_{red}=0.995$, $\phi =1$, and $N_{\mathrm{iter}}=600$. As for the SIRT-TV algorithm, ${\mathit{\lambda }}_{TV}=1.0$. ${\mathit{\lambda }}_{TV}$ is a positive factor of TV. As for the SIRT-WTDM algorithm, $\omega =0.00035$. As for the SIRT-MWTDM algorithm, $\gamma =0.6$ and $k=100$. The reconstructed images with different numbers of projections obtained by using the five reconstruction algorithms are shown in Figure 3.

From Figure 3, we can see that as the number of projections increases, the clarity and smoothness of the images reconstructed by the five algorithms are gradually improved, and the quality of reconstructed images is significantly improved when there is no noise in the projections. However, compared with the other four algorithms, there are obvious artifacts in the image reconstructed by using the FBP algorithm. The reconstructed images of SIRT, SIRT-TV, and SIRT-WTDM are smoother than that of FBP algorithm from 30 projections. The edges of the image reconstructed by using the SIRT-MWTDM algorithm are sharper than those of other algorithms. The boundaries of the three ellipses at the bottom of the image reconstructed by using SIRT-MWTDM algorithm are clearer than those of other algorithms. For 60 noise-free projections with 600 iterations, the image reconstructed by using SIRT-MWTDM method almost has no difference f the original image in visual effect. For 60 noisy projections, it can be seen that the image reconstructed by SIRT-MWTDM is less noisy. The intervals among the three ellipsoids are more obvious at the bottom of the reconstructed image.

SSIM, PSNR, and RMSE are calculated for the images reconstructed by using the five algorithms with different numbers of projections, and the results are shown in Figure 4. We first analyze the evaluation index of the image reconstructed by using noise-free projections. As can be seen from the SSIM histogram in Figure 4a, as the number of projections increases, the SSIM values of the image reconstructed by using each algorithm are improved. In the case of the same number of projections, the SIRT-MWTDM algorithm has the highest SSIM for reconstructed images. The SSIM value of the image reconstructed by the SIRT-MWTDM algorithm is 0.9878 for 30 projections, which also indicates that the proposed algorithm is superior to other algorithms in retaining structural information of the reconstructed image. As can be seen from the PSNR histogram in Figure 4b, PSNR and SSIM have similar characteristics. As the number of projections increases, the PSNR of the images reconstructed by using each algorithm is also improved. The PSNR value of the image reconstructed by the SIRT-MWTDM algorithm is the largest with the same number of projections. The PSNR value of the image reconstructed by using the SIRT-MWTDM algorithm is 34.3259 for 30 projections, which indicates that the proposed algorithm is superior to other algorithms in noise reduction. As can be seen from the RMSE histogram in Figure 4c, the variation trend of RMSE is opposite to that of SSIM and PSNR. The RMSE values of the image reconstructed by using SIRT-MWTDM are the smallest with the same number of projections compared with the other four algorithms. The RMSE values of each algorithm decrease with the increase in projection number. The RMSE value of the image reconstructed by the SIRT-MWTDM algorithm is 0.0145 for 30 projections, which indicates that the difference between the reconstructed image and the reference image is smaller and the image quality is higher. Then, we analyze the evaluation index of the image reconstructed by using the noise projection data. As can be seen from the SSIM and PSNR histogram, when the same number of noisy projections is used for reconstruction, both SSIM and PSNR values of the image reconstructed by using SIRT-MWTDM are the largest. As can be seen from the RMSE histogram, the RMSE value of the image reconstructed by using SIRT-MWTDM is the smallest when the same number of noisy projections is used for reconstruction. Therefore, from the reconstructed results obtained by using projection data with and without noise, it can be concluded that the objective evaluation indexes of images reconstructed by using our algorithm are the best compared with the other four algorithms under the same conditions.

To evaluate the performance of the SIRT-MWTDM algorithm, we plot the profile of the blue line position in Figure 2a as shown in Figure 5. From left to right, the number of projections is 30, 60, and 60+. The profiles demonstrate that the amplitude of the FBP and SIRT algorithms considerably fluctuates and differs from the true value of the pixel. SIRT-TV and SIRT-WTDM perform much better than the FBP and SIRT algorithms, but there is still a certain gap between their amplitude and the true value. However, the result of the SIRT-MWTDM algorithm matches that of the original model and is much closer to the true value. It also shows the effectiveness of SIRT-MWTDM algorithm in sparse-view CT reconstruction.

3.2.2. Circular Phantom

In the circular phantom simulation experiment, we set the parameters of the algorithms as follows: ${\lambda }_{\mathrm{SIRT}}=1.6$, ${\lambda }_{red}=0.995$, $\phi =1$, and $N_{\mathrm{iter}}=600$. As for the SIRT-TV algorithm, ${\lambda }_{TV}=1.0$. As for the SIRT-WTDM algorithm, $\omega =0.00085$. As for the SIRT-MWTDM algorithm, $\gamma =0.6$ and $k=90$. The reconstructed images obtained by using the five algorithms with different numbers of projections are shown in Figure 6.

From Figure 6, we can see that in the case of the same number of projections, the images reconstructed by using the SIRT-MWTDM algorithm have the best reconstruction effect. The outline of each circle from the reconstructed image can be seen clearly by using the SIRT-MWTDM algorithm with 30 projections. Even if there is Poisson noise in 60 projections, we can still see the position and contour of each circle. From the visual point of view, the reconstructed image obtained by our algorithm is smoother and has less noise for 60+ projections.

In order to quantitatively analyze the reconstructed results of the different algorithms, the histograms of three indexes (SSIM, PSNR, and RMSE) of the reconstructed image of the circular phantom are shown in Figure 7. From left to right, the histograms respectively represent the SSIM, PSNR, and RMSE of the reconstruction results of different algorithms for different numbers of projections. The SSIM and PSNR of five different algorithms show that the SIRT-MWTDM algorithm has the highest SSIM and PSNR in the same number of projections. The SIRT-MWTDM algorithm has the lowest RMSE in the same condition of projections. Therefore, Figure 7 shows that the image reconstructed by using the SIRT-MWTDM algorithm has the strongest noise suppression and the best performance.

To further investigate whether the performance of the proposed algorithm is better, we plot the profile of the red line position in Figure 2b as shown in Figure 8. From left to right, the number of sparse-view angles is 30, 60, and 60+. Here, 60+ indicates that Poisson noise is added to 60 projections. It can be seen from Figure 8 that the profile of the image reconstructed by using the SIRT-MWTDM algorithm is the closest to the original profile. In addition, the algorithm well protects the edge of the original image.

3.3. Neutron Experiment

The neutron experiments were carried out on a neutron imaging test station at end of the CNGB guide of the China Advanced Research Reactor (CARR) located at the China Institute of Atomic Energy. In order to simulate the neutron imaging nondestructive testing process of spent fuel, a 20 cm simulated nuclear fuel rod was used to obtain neutron projection data by using an indirect neutron CT device based on a neutron IP plate [33]. During the experiment, the operation power of CARR was 30 MW, and the neutron flux at the sample was 1 × 10⁸/cm²/s with a collimation ratio L/D of approximately 100. In the experiment, the 20 cm × 25 cm neutron IP plates were used to obtain the neutron projection images, and the imaging data were obtained by scanning the IP plates with the Typhoon FLA 7000 IP laser scanning imager. The neutron projection images of the nuclear fuel element were collected every 1° in the range of $[0^{\mathrm{o}},{180}^{\mathrm{o}}]$. Each neutron IP plate collected 5 projection data, and a total of 36 neutron IP plates with 180 projection data were collected. A normalized neutron projection image from one angle is shown in Figure 9a. The sinogram of an axial layer is shown in Figure 9b. The low quality of the normalized neutron projection image is due to the noise produced by gamma rays on the neutron IP plates, the inconsistency of the neutron IP plates in the experiment, and the registration error of the projection images in the angle segmentation. Since the nuclear fuel element is only close to the middle part in the projected image, the image is clipped before reconstruction to leave only the image of the nuclear fuel element. The size of the useful image is $5000\times 200$.

FBP algorithm is used to restore slice image from 180 projections, and the reconstructed image is used as a reference image to compare with the reconstructed results obtained by using other algorithms (FBP, SIRT, SIRT-TV, SIRT-WTDM, and SIRT-MWTDM) with 60 projections.

The reconstructed reference image is shown in Figure 10a, and it is normalized. The reconstructed results of other algorithms are also normalized. The reconstructed results of the five algorithms are shown in Figure 10b–f. The parameters of the iterative algorithm are as follows: ${\lambda }_{\mathrm{SIRT}}=0.8$, ${\lambda }_{red}=0.995$, $\phi =1$, and $N_{\mathrm{iter}}=100$. As for the SIRT-TV algorithm, ${\lambda }_{TV}=1.0$. As for the SIRT-WTDM algorithm, $\omega =0.00003$. As for the SIRT-MWTDM algorithm, $\gamma =0.3$ and $k=20$.

From the reconstructed results in Figure 10b, it is observed that severe artifacts exist in this image due to the limitation of the number of projections. The interior of the reconstructed image is rougher because of noise. The SIRT algorithm provides a better reconstructed image compared with FBP, but the streak artifacts are still present as shown in Figure 10c. As shown in Figure 10d,e, a large number of artifacts are suppressed by using the SIRT-TV and SIRT-WTDM algorithms. However, the effect of the SIRT-TV and SIRT-WTDM algorithms on noise suppression is not as good as that of the SIRT-MWTDM algorithm. Visually, the SIRT-MWTDM algorithm performs better in terms of noise suppression and artifact reduction compared with the other algorithms. Moreover, the SIRT-MWTDM algorithm is superior to other methods and produces an excellent-quality image in terms of maintaining detailed structure. It preserves most structural information without any physical deformation.

The RMSE, PSNR, and SSIM are used to carry out quantitative analysis on the reconstructed results obtained by using the five algorithms for 60 projections. The results are shown in Figure 11. It can be seen from Figure 11 that the RMSE value of the image reconstructed by SIRT-MWTDM algorithm is the smallest, which indicates the minimum error between the reconstructed image and the reference image. The PSNR value of the image reconstructed by using SIRT-MWTDM algorithm is the largest, which indicates that the reconstructed image is smoother and the algorithm has a stronger ability to suppress noise. The SSIM value of the reconstructed image obtained by SIRT-MWTDM algorithm is the largest, which indicates that the reconstructed image retains more structural information.

We also plot the profile in Figure 12, which shows the profile marked by the horizontal yellow line in Figure 10a. By comparison, the SIRT-MWTDM algorithm has advantages in that the profile is much closer to the reconstructed results obtained by using the FBP algorithm and 180 projections. In conclusion, real neutron reconstructed images indicate that our algorithm has good performance in artifact reduction and noise suppression.

AVIZO software was used to measure the area of the high-brightness substance indicated by the red arrow in Figure 10, and the measurement results are shown in

Table 1. The area of high-brightness substance.

Image	Figure 10a	Figure 10b	Figure 10c	Figure 10d	Figure 10e	Figure 10f
Area (mm2)	0.3850	0.2075	0.3300	0.3375	0.3450	0.3775

. As can be seen from Table 1, the area measured from Figure 10f is closest to that of the reference image Figure 10a.

The reconstruction of the nuclear fuel element was performed from all axial positions by using the proposed algorithm. A three-dimensional (3D) representation of several series of slices is shown in Figure 13. The AVIZO software was used to generate the 3D representation. To show the internal structure more clearly, the 3D image is displayed in color. Some of the pixels in the 3D image of the nuclear fuel element were dissected out to see the internal structure. The boundary of the nuclear fuel element and the vertical hole can be seen clearly in the diagram. The different colors of the image represent different types of substances. The red pixels represent uranium, which has a high neutron attenuation factor, and we can see the spatial distribution of uranium in the nuclear fuel element. A high-quality 3D image is very meaningful for researchers in the study of changes in the internal substances of nuclear fuel elements.

4. Discussion

The purpose of this paper is to obtain high-quality reconstructed images through incomplete real neutron projections. Because the traditional FBP algorithm has difficulty recovering the real structural information of the sample due to the high noise and a small number of neutron projection data, we consider using an iterative algorithm to solve this problem.

When the SIRT-MWTDM algorithm is used to reconstruct an image, some parameters need to be set, such as ${\lambda }_{\mathrm{SIRT}}$, $N_{\mathrm{iter}}$, $\gamma$, and $k$. In our research, we found that if the parameter ${\lambda }_{\mathrm{SIRT}}$ fluctuates in a small range, the quality of the reconstructed image does not change much. In other words, the SIRT-MWTDM algorithm is not very sensitive to this parameter. We can adjust ${\lambda }_{\mathrm{SIRT}}$ through experience and reconstructed effect. $N_{\mathrm{iter}}$ is usually set between 100 and 200 because as the number of iterations increases, the speed of image quality improvement will become slower and slower. Generally, after 100 or 200 iterations, a better reconstruction result of the algorithm can be achieved. However, more experience is required to set the values of $\gamma$ and $k$. According to the quality of the reconstructed results, we need to continuously adjust the values of these two parameters until a satisfactory result is obtained, which would be the shortcoming of the proposed algorithm in this paper. For the reconstruction of new samples, it may take some time to adjust the parameters $\gamma$ and $k$. The other algorithms used in the article use fewer parameters than the proposed algorithm. Therefore, it takes less time for other algorithms to adjust parameters than the proposed algorithm. We believe it is worthwhile to spend some time adjusting the parameters of the proposed algorithm to achieve high-quality reconstructed images.

We compare the convergence speed of the SIRT-MWTDM method with other algorithms. The RMSE and SSIM curves of the reconstruction process in Figure 3 for 30 Shepp–Logan projections are shown in Figure 14. The RMSE curve for each iteration is shown in Figure 14a. SIRT-MWTDM algorithm has the fastest descent speed at the beginning and the lowest final RMSE value from Figure 14a. The results show that the SIRT-MWTDM algorithm has a faster convergence speed and the gray value of the reconstructed image is closer to the ideal gray value. Figure 14b presents the SSIM curve for each iteration. SIRT-MWTDM algorithm has the fastest rising speed at the beginning and the largest final SSIM value from Figure 14b. When the three iterative algorithms (SIRT, SIRT-TV, SIRT-WTDM) are used to reconstruct the Shepp–Logan phantom from spare projections, they take about 100 iterations to achieve convergence, while the SIRT-MWTDM algorithm in this paper only needs about 50 iterations to achieve convergence.

The RMSE and SSIM curves of the reconstruction process in Figure 6 for 30 projections of the circular phantom are shown in Figure 15. Similarly, we can find from Figure 15 that as the number of iterations increases, the SIRT-MWTDM algorithm has the smallest RMSE and the largest SSIM. All algorithms reach convergence in about 50 iterations. Comparing Figure 14 and Figure 15, it can be found that the number of iterations to achieve convergence is not the same when the same algorithm is used to reconstruct different phantoms. In contrast, our method requires a small number of iterations to reach convergence. The image reconstructed by our method in this paper can obtain the smallest RMSE and the largest SSIM for the same number of projections.

5. Conclusions

MRP regularization can be used to suppress local non-monotonic structures and eliminate the influence of isolated noise. WTDM regularization not only considers the sparsity constraint of the image gradient, but also considers the continuity constraint of the image gradient, which can effectively protect the gradient information in different directions of the image. Therefore, we combine MRP and WTDM as the regularization constraint of the SIRT-MWTDM algorithm.

The quality of the reconstructed image is often affected by noise in the actual acquired data. To simulate the real image more realistically, Poisson noise was added to the simulated sinogram. Different numbers of projections were used to reconstruct images. We used subjective evaluation and three objective indicators (PSNR, SSIM, and RMSE) to evaluate and analyze the reconstructed results obtained by using five algorithms (FBP, SIRT, SIRT-TV, SIRT-WTDM, and SIRT-MWTDM). The SIRT-MWTDM algorithm shows advantages compared with the other algorithms both in subjective vision and objective evaluation. Moreover, the profiles from the interesting positions of the reconstructed image also show that the modified algorithm in this paper can better retain the details and suppress noise. The profile of the image reconstructed by the modified algorithm is closer to the reference image. Finally, it can be concluded that the SIRT-MWTDM algorithm shows great advantages both in simulated data and real neutron projection data when compared to the other algorithms in this paper.