TV algorithm

Rudin, Osher, and Fetamin proposed an ROF model using the L1 parameterization of the image gradient as regularization [24]. The model uses TV regularization to minimize the total variance and improve the image quality. The ROF model uses the L1 norm in the minimization process, which maintains the nonlinear characteristics of the image and thus allows more image edge structures to be retained. The expression for this model is as follows: $f_{\text{reg}}=\arg \mathrm{ min}\left\{{\Vert p-wf\Vert }_2^2+\lambda {\Vert f\Vert }_{\text{TV}}\right\},$ (1) where ${\Vert p-wf\Vert }_2^2$ is a data fidelity item and ${\Vert f\Vert }_{\text{TV}}$ is the regularization term, which is the TV norm of the image. The TV norm is the L1 norm of the image gradient, and can be expressed as follows: ${\Vert f\Vert }_{\text{TV}}={\displaystyle \sum \sqrt{{\left|{\nabla }_{x}f\right|}^2+{\left|{\nabla }_{y}f\right|}^2+{\left|{\nabla }_{z}f\right|}^2}},$ (2) where f denotes a 3D image, ▽f_i represents the gradient along the i direction, and |·| represents the complex modulus.

Split Bregman solved total variation (SBTV)

The Bregman distance [33], also known as the Bregman scatter, between two points u and v, is defined by the following equation: $D_E^p\left(u,v\right)=E\left(u\right)-E\left(v\right)-\langle p,u\text{-}v\rangle ,$ (3) where E(u) is a convex function, p denotes the subgradient of E(u), and <,> denote the inner product.

In the optimization problem, the Bregman distance is introduced as a regularization term. The Bregman iterative regularization method can be expressed using the imaging principle formula as follows: $u^{k+1}=\underset{u}{\min }{D}_E^{p^k}\left(u,u^k\right)+\frac{u}{2}{\Vert Au-f\Vert }_2^2,$ (4) since Eq. (4) is a convex optimization problem. Substituting this into Eq. (5), the following results can be obtained using the optimal conditions: $p^{k+1}=p^k-\mu A^{\text{T}}\left(A{u}^{k+1}-f\right).$ (5)

The Bregman iterative method comprises Eq. (4) and Eq. (5). Eq. (4) optimizes the subproblem, and Eq. (5) updates the search direction of the iterative solution. However, owing to the difficulty of the subgradient calculation, the two equations above are only warped according to the nature of recursion. Therefore, Bregman’s iterative method is obtained in a simplified form [34] as $\begin{array}{l}{u}^{k+1}=\underset{u}{\min }E\left(u\right)+\frac{u}{2}{\Vert Au-f-b\Vert }_2^2\\ b^{k+1}=b^k+\left(A{u}^{k+1}-f\right).\end{array}$ (6)

When used to reconstruct NCT images, the reconstruction algorithm did not provide satisfactory reconstruction results. Therefore, the image must be preprocessed, and generally good results are achieved using the denoising TV method. The objective function of image denoising using the split-Bregman method is expressed as follows: $\underset{f}{\min }\left|{\nabla }_{x}f\right|+\left|{\nabla }_{y}f\right|+\left|{\nabla }_{z}f\right|+\frac{\mu }{2}{\Vert f-f^{it}\Vert }_2^2.$ (7)

Solving Eq. (7) using the split Bregman method transforms the problem into three easily solvable subproblems. Therefore, transforming Eq. (7) into a constrained problem, it can be expressed as: $\underset{f,d_x,d_y,d_z}{\min }\left|d_x\right|+\left|d_y\right|+\left|d_z\right|+\frac{\mu }{2}{\Vert f-f^{it}\Vert }_2^2,$ (8) where d_x=▽_xf, d_y=▽_yf, d_z=▽_zf Subsequently, using Bregman iteration, we solved Eq. (7), which can be further expressed as $\begin{array}{l}\underset{f,d_x,d_y,d_z}{\min }\left|d_x\right|+\left|d_y\right|+\left|d_z\right|+\frac{\mu }{2}{\Vert f-f^{it}\Vert }_2^2\text{+}\frac{\lambda }{2}{\Vert d_x^k-{\nabla }_{x}d-b_x^k\Vert }_2^2\\ +\frac{\lambda }{2}{\Vert d_y^k-{\nabla }_{y}d-b_y^k\Vert }_2^2+\frac{\lambda }{2}{\Vert d_z^k-{\nabla }_{z}d-b_z^k\Vert }_2^2\\ f^{k+1}=\underset{f}{\min }\frac{\mu }{2}{\Vert f-f^{it}\Vert }_2^2\text{+}\frac{\lambda }{2}{\Vert d_x^k-{\nabla }_{x}f-b_x^k\Vert }_2^2\\ +\frac{\lambda }{2}{\Vert d_y^k-{\nabla }_{y}f-b_y^k\Vert }_2^2+\frac{\lambda }{2}{\Vert d_z^k-{\nabla }_{z}f-b_z^k\Vert }_2^2\\ d_i^{k+1}=\underset{d_i}{\min }{\Vert d_i\Vert }_1+{\Vert d_i-{\nabla }_i{f}^{k+1}-b_i^k\Vert }_2^2.\end{array}$ (9)

As we have observed, the split Bregman approach differentiates the problem by separating the L1 part. We applied the Gauss-Seidel method [35] to solve f in a simpler manner, as follows: $\begin{array}{c}{f}_{g,j,v}^k=\frac{\lambda }{\mu +6\lambda }(f_{g+1,j,v}^k+f_{g-1,j,v}^k+f_{g,j+1,v}^k\\ +f_{g,j-1,v}^k+f_{g,j,v+1}^k+f_{g,j,v-1}^k\\ +d_{x,g-1,j,v}^k-d_{x,g,j,v}^k+d_{y,g,j-1,v}^k\\ -d_{y,g,j,v}^k+d_{z,g,j,v}^k-d_{z,g,j,v}^k\\ -b_{x,g-1,j,v}^k+b_{x,g,j,v}^k-b_{y,g,j-1,v}^k+b_{y,g,j,v}^k\\ -b_{z,g,j,v-1}^k+b_{z,g,j,v}^k)+\frac{\lambda }{\mu +6\lambda }{f}_{g,j,v}.\end{array}$ (10)

Additionally, since the parts of $d$. e unrelated, they can be separated into dx, dy, and dz. Thus, we employ the shrinkage operator to obtain the optimal solution [36], which can be written as follows: $d_i^{k+1}=shink\left({\nabla }_i{f}^{k+1}+b_i^k,\frac{1}{\lambda }\right)=\frac{{\nabla }_i{f}^{k+1}+b_i^k}{\left|{\nabla }_i{f}^{k+1}+b_i^k\right|}\times \max \left(\left|{\nabla }_i{f}^{k+1}+b_i^k\right|-\frac{1}{\lambda },0\right).$ (11)

3D Anisotropic diffusion filtering algorithms (3D-ADF)

Sun et al. extended two-dimensional anisotropic diffusion filtering to three dimensions and proposed a three-dimensional anisotropic diffusion filtering algorithm that can retain image edge structure information, as well as perform denoising [37]. $\frac{\partial f\left(i,l,r;t\right)}{\partial t}=div\left[q\left(i,l,r;t\right)\nabla f\left(i,l,r;t\right)\right],$ (12) where f denotes the three-dimensional image, ∇ represents the gradient operator, and div denotes the diffusion coefficient. Current research on anisotropic diffusion-filtering algorithms focuses on the selection of c(x).

In the algorithm proposed by Tang, c(x) is expressed as follows: $q\left(i,l,r;t\right)=\frac{1}{1+\left[\frac{{\left|\nabla f\left(i,l,r;t\right)\right|}^2}{f{\left(i,l,r;t\right)}^2}-p^2\right]/\left(p^2\left(1+p^2\right)\right)}.$ (13)

The equation for the diffusion coefficient of 3D anisotropic diffusion is $q\left(h\right)=\frac{1}{1+\left[h^2\left(i,l,r;t\right)-h_0^2\left(t\right)\right]/\left[h_0^2\left(t\right)\left(1+h_0^2\left(t\right)\right)\right]}$ (14) $h_0\left(t\right)=\frac{\sqrt{\mathrm{var}\left[z\left(t\right)\right]}}{\overline{z\left(t\right)}},$ (15) where var[z(t)] denotes the intensity variance, and $\overline{z\left(t\right)}$ represents the average value of the speckle area at position t. The expression for the instantaneous coefficient of variation h(i,l,r;t) in three dimensions is as follows: $h\left(i,l,r;t\right)=\sqrt{\frac{\left(1/3\right){\left(\left|\nabla f\right|/f\right)}^2-{\left(1/6\right)}^2{\left({\nabla }^{2}f/f\right)}^2}{{\left(1+\left(1/6\right)\left({\nabla }^{2}f/f\right)\right)}^2}}.$ (16)

These algorithms not only remove noise, but also retain more information about the image edge structure. $\begin{array}{c}\frac{\partial f\left(i,l,r;t\right)}{\partial t}=div\left[q\left(i,l,r;t\right)\nabla f\left(i,l,r;t\right)\right]\\ f\left(i,l,r;0\right)=f_0\left(i,l,r\right).\end{array}$ (17)

In order to solve Eq. (17), it was discretized, and the discrete expression of the 3D anisotropic diffusion is obtained as in Eq. (18): $\begin{array}{l}{f}_{g,j,k}^k=f(gh,jh,vh,k\text{Δ}t)\\ \{\begin{array}{l}{f}_{g,j,v}^{k+1}=f_{g,j,v}^k+\frac{\text{Δ}t}{6}{d}_{g,j,v}^k\\ d_{g,j,v}^k=\frac{1}{h^2}\left[\begin{array}{l}{q}_{g+1,j,v}^k\left(f_{g+1,j,v}^k-f_{g,j,v}^k\right)+q_{g,j,v}^k\left(f_{g-1,j,v}^k-f_{g,j,v}^k\right)\\ +q_{g,j+1,v}^k\left(f_{g,j+1,v}^k-f_{g,j,v}^k\right)\end{array}\right].\\ +\frac{1}{h^2}\left[\begin{array}{l}{q}_{g,j,v}^k\left(f_{g-1,j,v}^k-f_{g,j,v}^k\right)+q_{g,j,v+1}^k\left(f_{g,j,v+1}^k-f_{g,j,v}^k\right)\\ +q_{g,j,v}^k\left(f_{g,j,v-1}^k-f_{g,j,v}^k\right)\end{array}\right]\end{array}\end{array}$ (18)

OS-SART algorithm

The OS-SART algorithm is a variant of the SART algorithm. By repeatedly breaking up the projection data into ordered subgroups for each subset, the computing efficiency of SART is increased. At each iteration, the algorithm reconstructs a subset of the image based on the corresponding projected subset. This procedure was repeated for each subgroup, and the results were combined to create a final reconstructed image. By using ordered subsets, OS-SART creates a tradeoff between the computation time and image quality. It provides faster reconstruction speeds than conventional SART, while still achieving better image quality. By varying the number of subsets and iterations, the equilibrium between the reconstruction speed and image accuracy can be regulated. The formula is given by the equation below: $f_j^{\left(k+1\right)}=f_j^k+\frac{{\displaystyle \sum {}_{i\in S_i}}\left\{{\lambda }_k\frac{p_i-{\displaystyle {\sum }_{m=0}^{\text{M-1}}{w}_{im}{f}_{im}^{\left(k\right)}{w}_{ij}}}{{\displaystyle {\sum }_{m=0}^{\text{M-1}}{w}_{im}}}\right\}}{{\displaystyle \sum {}_{i\in S_i}{w}_{ij}}},$ (19) where i represents the i-th ray at a particular projection angle, j refers to the pixel value in the i-th ray, k represents the number of iterations, and λ is a relaxation factor.

Proposed algorithm

Based on the OS-SART image reconstruction and 3D ADF image processing algorithms, we presented a sparse reconstruction algorithm (OS-SART-ADF) for NCT. The algorithm comprises two loops: the outer loop employs OS-SART for image reconstruction, and the inner loop employs ADF for image denoising.

Quantitative Evaluation Index

The correlation coefficient (CC) is a commonly used metric in image evaluation that utilizes the covariance between the pixel values of two images to measure the correlation between them using the following formula: $\text{CC=}\frac{{\displaystyle {\sum }_{k=1}^M\left(\varpi \left(y\right)-\overline{\varpi }\right)\left({\varpi }_{\text{true}}\left(y\right)-{\overline{\varpi }}_{\text{true}}\right)}}{{\left\{{\displaystyle {\sum }_{k=1}^M{\left(\varpi \left(y\right)-\overline{\varpi }\right)}^2}{\displaystyle {\sum }_{y=1}^M{\left({\varpi }_{\text{true}}\left(y\right)-{\overline{\varpi }}_{\text{true}}\right)}^2}\right\}}^{1/2}},$ (20) where ${\overline{\varpi }}_{\text{true}}$ represents the measurement image, M denotes the voxel number, and ${\overline{\varpi }}_{\text{true}}$ represents the reconstructed value at voxel y.

The mean structural similarity (MSSIM) was used to measure the structural similarity between two images. This takes into account information such as contrast, brightness, and structure of the image to provide a more comprehensive evaluation.

Brightness information: $l\left(\epsilon ,\eta \right)=\frac{2u_{\epsilon }{u}_{\eta }+C_1}{u_{\epsilon }^2+u_{\eta }^2+C_1}$

Contrast information: $c\left(\epsilon ,\eta \right)=\frac{2{\sigma }_{\epsilon }{\sigma }_{\eta }+C_2}{{\sigma }_{\epsilon }^2+{\sigma }_{\eta }^2+C_2}$

Structure Information: $s\left(\epsilon ,\eta \right)=\frac{{\zeta }_{\epsilon \eta }+C_3}{{\zeta }_{\epsilon }+{\zeta }_{\eta }+C_2}$

where ε,η represents the reference image and measurement image, separately; u_ε and u_η are the averages of the image, σ_ε and σ_η are the standard deviations of the image, ζ_εη is the covariance of the image, and C₁, C₂, and C₃ are constants. Therefore, the MSSIM formula is as follows: $\text{MSSIM=}\frac{1}{\text{M}}{\displaystyle \sum _{t=1}^{\text{M}}{\left[l\left({\epsilon }_t,{\eta }_t\right)\right]}^{\alpha }}\times {\left[c\left({\epsilon }_t,{\eta }_t\right)\right]}^{\beta }\times {\left[s\left({\epsilon }_t,{\eta }_t\right)\right]}^{\gamma }.$ (21)

The root mean square error (RMSE) is used to measure the difference between two images; it is a measure of similarity between two images obtained by calculating the difference between their corresponding pixel values. $\text{RMSE=}\sqrt{\frac{{\displaystyle {\sum }_{n=1}^{\text{N}}\left(\overline{f_n}-f_n\right)}}{\text{N}}},$ (22)

Here, $\overline{f_n}$ denotes a voxel in the measurement image, f_n denotes a voxel in the sample, and N denotes the number of voxels.

The universal quality index (UQI) is an objective metric used to measure the quality of an image, and is calculated by comparing the differences in structure, brightness, and contrast between a reference image and measurement image. $\text{UQI}=\frac{2\mathrm{cov}\left(\overline{\mu },\mu \right)}{{\overline{\rho }}^2+{\rho }^2}\cdot \frac{2\overline{\mu }\mu }{{\overline{\mu }}^2+{\mu }^2},$ (23)

In the above equation, cov denotes the pixel value covariance of the reference and measurement images, $\overline{\mu }$ and μ denote the mean pixel values of the reference and measurement images, respectively, and ${\overline{\rho }}^2\text{ and }{\rho }^2$ represent the standard deviation of the pixel values of the reference and measurement images, respectively.

Neutron noise model

Photon statistical noise is a statistical fluctuation in the number of photons received per pixel by the detector owing to the limited number of rays. The photon statistical noise obeys a Poisson distribution and manifests itself as random variations in the image, particularly in low-dose scans. Electronic noise is generated by the detectors and electronics in the CT system and can affect the quality of the image. The electronic noise typically follows a Gaussian distribution. Additionally, a variety of noise is present in PET, MRI, and other imaging systems, and various algorithms and filters must be used in the reconstruction of images from these imaging systems to minimize the effects of noise and retain useful information [38-39]. The following equation can be used in NCT imaging to represent the noise model of the system: $\overline{L}=Poisson(L)+Gaussian\left(m_{\text{ic}}+{\delta }_{\text{ic}}^2\right),$ (24) where $\overline{L}$ denotes noisy transmission data, $L$ is the mean number of photons, $m_{\text{ic}}$ and ${\delta }_{\text{ic}}^2$ are the mean and variance of the noise of the electronic device, respectively. The value $m_{\text{ic}}$ is generally set to zero when commissioning the system. The variance ${\delta }_{\text{ic}}^2$ is assessed by dark current measurements [40]. In this study, three parameters, the photon incident flux 1×10⁵, $m_{\text{ic}}=0$, and ${\delta }_{\text{ic}}^2=0.25$ were set based on the characteristics of the NCT system.

Digital 3D Shepp-Logan Model Experiment

By reconstructing the digital 3D Shepp-Logan model with dimensions of 256×256×256, we compare and assess the performance variances of four algorithms: filtered back-projection (FBP), OS-SART-TV, OS-SART-SBTV, and OS-SART-ADF, by reconstructing the digital 3D Shepp-Logan model of size 256×256×256. A slice image of the 3D Shepp-Logan mode at z=128 is shown in Fig. 1. The three iterative algorithm parameters are specified as follows: OS-SART-TV, ${\lambda }_{\text{OS-SART}}=2.2$, ${\text{Niter}}_{\text{OS-SART}}=50$, ${\lambda }_{\text{TV}}=200$, AND $\text{Niter=300}$. As for the OS-SART-SBTV algorithm, ${\lambda }_{\text{OS-SART}}=2.2$, ${\text{Niter}}_{\text{OS-SART}}=50$, ${\lambda }_{\text{TV}}=0.0032$, ${\text{Niter}}_{\text{TV}}=150$. As for the OS-SART-ADF algorithm (Table 1), ${\lambda }_{\text{OS-SART}}=2.2$, ${\text{Niter}}_{\text{OS-SART}}=50$, ${\lambda }_{\text{ADF}}=0.002$, ${\text{Niter}}_{\text{ADF}}=300$.

Fig. 1Full-screen View

3D Shepp-Logan slice image at z=128

The images produced by the four algorithms for reconstructing the different projection angles are shown in Fig. 2. The findings of all the algorithms show that an increase in the projected angle count significantly improves the smoothness of the reconstructed image and sharpness of the edge structure. When comparing the iterative reconstruction results, the reconstructed image of the FBP algorithm revealed significant artifacts caused by the algorithm’s design defects. By comparison, the images reconstructed using the three iterative algorithms had clearer edge-structure information. From the region of interest in Fig. 2, it is clear that although OS-SART-TV can effectively reduce artifacts, the fine structure of the image is lost more significantly, resulting in the reconstructed image being too smooth. When the number of sparse views is 25, OS-SART-SBTV can also reconstruct clearer structural information of the image, but there is still ambiguity, in contrast to OS-SART-ADF. By comparing the reconstructed images from the different algorithms, it was observed that the OS-SART-ADF algorithm effectively eliminated image noise and clearly reconstructed the edge structure of the image.

Fig. 2Full-screen View

Sliced reconstructed images at z=128 for the four algorithms

We compared and analyzed the error images reconstructed using the five algorithms. The difference in pixel values between the original and reconstructed images at the same locations is reflected in the error image; the smaller the difference, the more structural information is included in the reconstructed image. From the error images in Fig. 3, the error image of FBP has the largest difference, indicating that the algorithm reconstructs the image with a large amount of structural feature information lost. In comparison with the analytical FBP algorithm, all these iterative algorithms performed well in removing artifacts and reconstructing the detailed structural information of the image. However, OS-SART-ADF had the least error in the image pixel values, indicating that it reconstructed the most detailed image.

Fig. 3Full-screen View

The slice of the error images at z=128

The reconstructed image contours of all algorithms are shown in Fig. 4 to provide an in-depth analysis of their reconstruction performance. From left to right, the projected views are 25, 35, and 35+, where 35+ denotes the image reconstructed with the added neutron noise model. The profile lines of the FBP-rebuilt image differ significantly from those of the reconstructed image, as shown in Fig. 4. Although OS-SART-TV and OS-SART-SBTV are superior to FBP, their reconstructed image profiles deviate slightly from the profiles of the reference image. In comparison, OS-SART-ADF reconstructs images whose profiles are most similar to the reference image. These results further demonstrate that OS-SART-ADF is highly appropriate for NCT.

Fig. 4Full-screen View

Profiles of reconstructed images using different algorithms:(a) horizontal profile and (b) vertical profile

The image quality of each algorithm was quantitatively evaluated using four image assessment indices, as shown in Fig. 5. We first analyzed all algorithms that reconstructed the image obtained with noiseless projection data. The image evaluation indices gradually decreased with the reduction in the projection angle according to the three image evaluation metrics CC, MSSIM, and UQI; the change trends of these metrics were identical. From the RMSE metric in Fig. 5, it can be seen that the OS-SART-ADF-reconstructed images had the minimum RMSE values under the same conditions. For instance, the RMSE of the OS-SART-ADF-reconstructed image is 0.0292 when the reconstructed projection view is 25, which demonstrates that the algorithm reconstructs the image with pixel values closest to the original values, and that the reconstructed image has the best quality. From the UQI metric in Fig. 5d, it can be seen that the UQI of all algorithms progressively increases with the number of projection views, which is opposite to the trend of the RMSE. OS-SART-ADF has the highest UQI for reconstructing images under identical conditions compared to all other reconstruction methods, indicating its superior performance. When OS-SART-ADF reconstructed the 25 projection views, the UQI, MSSIM, and CC of the reconstructed image were 0.9877, 0.9878, and 0.9887, respectively. These evaluation metrics of the OS-SART-ADF-reconstructed image were the highest among the three iterative algorithms. Subsequently, we analyzed the reconstructed images containing noisy projection data. It can be observed from the CC, MSSIM, and UQI in Fig. 5 that these metrics are maximized for the OS-SART-ADF algorithm-reconstructed images when the reconstructed projection views are the same, with values of 0.991, 0.9903, and 0.9901, respectively. The RMSE of the OS-SART-ADF-reconstructed images was minimized according to the RMSE in Fig. 5. As a result, when the reconstructed projection angles are the same, OS-SART-ADF has the best quality of reconstructed images, based on an examination of the reconstructed images.

Fig. 5Full-screen View

(Color online) Metrics of 3D Shepp-Logan model reconstructed images

3D Industrial Phantom

To further verify the performance of the OS-SART-ADF in real-world situations, we compared the reconstructed images of a 3D industrial model with dimensions of 256 × 256 × 256. A sliced image of the 3D industrial model at z=128 is shown in Fig. 6. By reconstructing the 3D industrial model, we evaluated the ability of each algorithm to retain structural information at the edges of the image.

Fig. 6Full-screen View

The 3D industrial phantom slice image at z=128

Fig. 7 shows the reconstructed images for each algorithm when reconstructing the 3D industrial model at different projection angles. Depending on the reconstructed results of the algorithm, the larger the projection angle, the clearer the edge structure of the reconstructed image. In contrast to the other algorithms, the reconstructed image of the FBP algorithm is blurred, and the information of the image edge structure is significantly lost. By analyzing the reconstruction results of each algorithm, the OS-SART-TV algorithm reconstructed a more delicate image structure that was excessively smooth. Among these iterative algorithms, OS-SART-ADF reconstructs the image with the clearest edge structure and least noise effect.

Fig. 7Full-screen View

Sliced images of the 3D industrial phantom at z=128 for the four algorithms

As shown in Fig. 8, we performed a quantitative analysis of the rebuilt results using the three iterative algorithms. Initially, the reconstructed images for each algorithm were analyzed using noiseless projection data. Fig. 8 shows that the CC, MSSIM, and UQI metrics of all the algorithms progressively increased as the projection angle increased. The three indicators exhibited similar trends. Based on the RMSE metric presented in Fig. 8, the OS-SART-ADF-reconstructed image had the lowest RMSE. For instance, when 25 projection views were reconstructed, the RMSE of OS-SART-ADF was 0.0079, indicating that the reconstructed image from this algorithm was the closest to the real image. As shown in Fig. 8d, the UQI of all algorithms progressively increases with an increase in the number of projected views, which is opposite to the trend of the RMSE. Among all reconstruction algorithms, the OS-SART-ADF-reconstructed images have the highest UQI, which indicates superior performance. For example, the UQI of OS-SART-ADF was 0.9978 when the number of reconstruction angles was 25. Subsequently, we performed a quantitative analysis of the reconstructed images containing noise. The CC, MSSIM, and UQI metrics in Fig. 8 show that the OS-SART-ADF algorithm reconstructed images with the highest metrics when reconstructing the same amount of projection data containing noise. Furthermore, the RMSE metric in Fig. 8 indicates that the reconstructed images of the OS-SART-ADF algorithm had the smallest RMSE value. Thus, a quantitative and intuitive analysis of each algorithm indicates that OS-SART-ADF reconstructs the best images among all tested algorithms with the same number of reconstruction.

Fig. 8Full-screen View

(Color online) Evaluation metrics for the 3D industrial phantom reconstructed images

Real neutron projection experiment

We compared the reconstruction capabilities of all algorithms using a clock model (Fig. 9) with projection data provided by Schillinger [41]. Figure 10 shows a slice image of the clock model at z=201. The following three-step process was applied to the neutron photographs to obtain valid neutron projection data: First, for the dark- and bright-field images, white noise in the image was removed using a 5×5 median filter. For the neutron photographic image, 3×3 median filtering was applied, which not only protected the detailed information of the neutron photographic image, but also removed white noise from the image. Second, the pixel values of the dark-field and bright-field images provided are averaged separately, and the averaged dark-field background image data are subtracted from the averaged bright-field image data so that the corrected bright-field image data can be obtained. Finally, the entire neutron photographic image data were subtracted from the averaged dark-field image data, normalized, and negatively notated. The non-negative data were then corrected to obtain the neutron projection data for the applied image reconstruction.

Fig. 9Full-screen View

(Color online) Clock model

Fig. 10Full-screen View

The clock model slices the image at z=201

As shown in Fig. 11, the image recreated using the FBP algorithm contained significant artifacts. The three iterative reconstruction algorithms eliminated artifacts better and reconstructed more image information than the FBP algorithm. Compared with the FBP algorithm, the reconstructed images from OS-SART-TV have fewer artifacts; however, because of the excessively smooth image edges of OS-SART-TV, a significant amount of the structural information of the images is lost. While the reconstructed image of OS-SART-SBTV is sharper than that of OS-SART-TV, it is still fuzzy compared to that from OS-SART-ADF. Therefore, OS-SART-ADF exhibits good performance in artifact suppression and noise reduction based on the visual observation of the reconstructed images.

Fig. 11Full-screen View

Slice image of the reconstructed image of the four algorithms at z=201

As shown in Fig. 12, we quantitatively analyzed the clock model images reconstructed using the three iterative algorithms. From the RMSE metric shown in Fig. 12, it can be seen that OS-SART-ADF has the minimum RMSE, which means that the reconstructed image of OS-SART-ADF is the closest to the reference image. Moreover, the CC, MSSIM, and UQI of OS-SART-ADF were the largest among all algorithms. In conclusion, the OS-SART-ADF applied to neutron CT 3D reconstruction performed excellently in fine-structure reconstruction, noise removal, and artifact suppression, which shows that the OS-SART-ADF algorithm is appropriate for neutron CT 3D reconstruction of sparse views.

Fig. 12Full-screen View

(Color online) Metrics for quantitative analysis of reconstructed clock model images