Low-Rank Matrix Approximation (LRMA) Algorithm

The LRMA algorithm is a method for approximating the representation of a high-dimensional matrix as a low-rank matrix [30]. In image processing, this algorithm simplifies the data by preserving their main features and structure.

We first analyzed the problem of estimating a low-rank matrix $X$ from noisy projection data. $Y=X+W,Y,X,W\in {\mathcal{R}}^{m\times n},$ (1) where W denotes the noise model. The LRMA problem can be defined using the following equation: $\arg \underset{x}{\min }\left\{\Psi \left(x\right)=\frac{1}{2}{\Vert Y-X\Vert }_F^2+\lambda {\displaystyle \sum _{i=1}^k\varphi }\left({\sigma }_i\left(X\right);a\right)\right\},$ (2) where k=min(m, n) and ${\sigma }_i\left(X\right)$ are the ith singular values of the matrix X, $\varphi$ is a regularizer that induces sparsity and may be non-convex.

Total Variation (TV) Algorithm

Iterative algorithms have an advantage over analytical algorithms under the condition of limited projection views; however, when severely under-sampled, the reconstructed image still has severe artifacts, even with iterative algorithms. Because the CS theory was proposed, CS-based approaches have been successfully applied to eliminate artifacts and other aspects. CS theory can only be applied to sparse images or images that can be represented sparsely. Generally, neutron tomography images are not sparse but can be sparsely expressed by a sparse transformation using the following definition equation: $\left|\stackrel{\rightharpoonup }{\nabla }{q}_{u,v,z}\right|=\sqrt{{\left(q_{u,v,z}-q_{u-1,v,z}\right)}^2+{\left(q_{u,v,z}-q_{u,v-1,z}\right)}^2+{\left(q_{u,v,z}-q_{u,v,z-1}\right)}^2},$ (3) where $q_{u,v,z}$ denotes the pixel value at position $\left(u,v,z\right)$.

The TV is the total pixel value in the image. $q_{\text{TV}}={\displaystyle \sum _{u,v,z}\left|\stackrel{\rightharpoonup }{\nabla }{q}_{u,v,z}\right|} ={\displaystyle \sum _{u,v,z}\sqrt{{\left(q_{u,v,z}-q_{u-1,v,z}\right)}^2+{\left(q_{u,v,z}-q_{u,v-1,z}\right)}^2+{\left(q_{u,v,z}-q_{u,v,z-1}\right)}^2}},$ (4)

Mathematically, TV-based image reconstruction can be expressed as follows: $\underset{q}{\min }{q}_{\text{TV}},s.t.\omega q=p$ (5)

Fast Gradient Projection (FGP) Algorithm

Beck and Teboulle introduced the FGP algorithm to solve the TV, which is derived as follows [31]: $\underset{x\in C}{\min }{\left|x-b\right|}_F^2+2\lambda \text{TV}\left(x\right).$ (6)

The non-smooth nature of TV results in the inability to directly solve Eq. (6). To solve this problem, Chambolle used a gradient-based algorithm for the solution process. Based on this algorithm, a constrained problem pair is constructed. Certain symbols are assumed. The $\psi$ represents the set of matrix pairs (f, q) that satisfy the following relations: $\begin{array}{l}{\left(f_{i,j}\right)}^2+{\left(q_{i,j}\right)}^2\le 1,\left(0\le i\le h,0\le j\le g\right)\\ \left|f_{i,g+1}\right|\le 1,0\le i\le \text{h}\\ \left|q_{h+1,j}\right|\le 1,0\le j\le g\end{array}$

$\kappa :{ℝ}^{\left(h-1\right)\times g}\times {ℝ}^{h\times \left(g-1\right)}\to {ℝ}^{h\times g}$ is linear operation. Its definition formula is as below. $\begin{array}{l}\kappa {\left(f,q\right)}_{i,j}=\left(f_{i,j}-f_{i-1,j}\right)+\left(q_{i,j}-q_{i,j-1}\right)\\ 1\le i\le h,1\le j\le g\end{array}$

The operator ${\kappa }^T$ $\kappa :{ℝ}^{h\times g}\to {ℝ}^{\left(h-1\right)\times g}\times {ℝ}^{h\times \left(g-1\right)}$ adjacent to $\kappa$ uses the following computational formula. ${\kappa }^T\left(o\right)=(f,q)$ where $f,q\in {ℝ}^{m\times n}$ ${ℝ}^{m\times n}$ are the matrices that defines the equation. $\begin{array}{l}{f}_{i,j}=o_{i,j}-o_{i+1,j},i=1,\dots ,h-1,j=1,\dots ,g\\ q_{i,j}=o_{i,j}-o_{i,j+1},i=1,\dots ,h,j=1,\dots ,g-1\end{array}$

P_C is an orthogonal projection operator of C.Thus, when C=B_l,u, the $P_{B_{l,u}}$ is given by $P_{B_{l,u}}{\left(O\right)}_{i,j}=\{\begin{array}{ll}l\hfill & o_{ij}<l\hfill \\ o_{ij}\hfill & l\le o_{ij}\le u\hfill \\ u\hfill & o_{ij}>u\hfill \end{array}$

Assuming this notation, we obtained the dual problem in Eq. (6). This approach can explain the connection between the dual and primal optimal solutions. To explain this relationship further, we introduced the following proposition. We assumed that (f, q)∈Ψ is the optimal solution to this problem. $\underset{\left(f,q\right)\in \psi }{\min }\left\{h(f,q)\equiv -{\Vert H_C\left(b-\lambda \kappa (f,q)\right)\Vert }_F^2+{\Vert b-\lambda \kappa (f,q)\Vert }_F^2\right\}.$ (7)

The expression for H_C(o) is given below. $H_C\left(o\right)=o-P_C\left(o\right),o\in {ℝ}^{h\times g}.$ (8)

h in Eq. (7) is a continuous differentiable function whose gradient is defined as follows: $\nabla h(f,q)=-2\lambda {\kappa }^T{P}_C\left(b-\lambda \kappa \left(f,q\right)\right),$ (9) $T(o,f,q)=Tr\left(\kappa {\left(f,q\right)}^{T}o\right).$ (10)

Consequently, the problem in Eq. (6) can be transformed into the following equation: $\underset{o\in C}{\min }\underset{\left(f,q\right)\in \psi }{\max }\left\{{\Vert o-b\Vert }_F^2+2\lambda Tr(\kappa {\left(f,q\right)}^{T}o\right\}.$ (11)

Because the objective function is concave in f, q and convex in o, it is possible to swap the maximum and minimum orders. Thus, we obtain the following formula: $\underset{\left(f,q\right)\in \psi }{\max }\underset{o\in C}{\min }\left\{{\Vert o-b\Vert }_F^2+2\lambda Tr(\kappa {\left(f,q\right)}^{T}o\right\}.$ (12)

It is able to be reformulated as below. $\underset{(f,q)\in \psi }{\max }\underset{o\in C}{\min }\left\{{\Vert o-(b-\lambda \kappa (f,q))\Vert }_F^2-{\Vert b-\lambda \kappa (f,q)\Vert }_F^2+{\Vert b\Vert }_F^2\right\}.$ (13)

Therefore, the optimal solution of the above equation can be obtained. $o=P_C\left(b-\lambda \kappa (f,q)\right)$ (14)

We can then solve the following dyadic problem by applying Eq. (14) to Eq. (13) and omitting the permanent term in Eq. $\underset{(f,q)\in \psi }{\max }\left\{{\Vert P_C\left(b-\lambda \kappa (f,q)\right)-(b-\lambda \kappa (f,q))\Vert }_F^2-{\Vert b-\lambda \kappa \left(f,q\right)\Vert }_F^2\right\}.$ (15)

Because our objective is to resolve the dual problem of Eq. (6), and the gradient is represented by Eq. (11), we introduced a gradient projection algorithm, which is commonly used for denoising problems. Because the norm of $oϵ{ℝ}^{h\times g}$ is Frobenius norm. For $(\text{f}\text{, }\text{q})ϵ{ℝ}^{\left(h-1\right)\times g}\times {ℝ}^{h\times \left(g-1\right)}$, it is expressed as the following equation. $\Vert (f,q)\Vert =\sqrt{{\Vert f\Vert }_F^2+{\Vert q\Vert }_F^2}$ (16) $Lh\le {\lambda }^2$ (17)

The projection in the set $\psi$ could be calculated simply. For $(f,q)$, the projection $P_{\psi }\left(f,q\right)$ can be given by $P_{\psi }\left(f,q\right)=\left(l^1,l^2\right)$. The expression for $lϵR^{h\times g}$ is given below. $\begin{array}{l}{l}_{i,j}^1=\{\begin{array}{c}\frac{f_{i,j}}{\text{max}\left\{1,\sqrt{{\left(f_{i,j}\right)}^2+{\left(q_{i,j}\right)}^2}\right\}}\\ \frac{f_{i,j}}{\text{max}\left\{1,\left|f_{i,j}\right|+\left|q_{i,j}\right|\right\}}\end{array}\\ l_{i,j}^2=\{\begin{array}{c}\frac{q_{i,j}}{\text{max}\left\{1,\sqrt{{\left(f_{i,j}\right)}^2+{\left(q_{i,j}\right)}^2}\right\}}\\ \frac{p_{i,j}^2}{\text{max}\left\{1,\left|f_{i,j}\right|+\left|q_{i,j}\right|\right\}}\end{array}\end{array}$ (18)

The image-denoising algorithm based on the FGP algorithm can be obtained by substituting the maximum values of the gradient equation, objective function, and Lipschitz constant into the gradient projection algorithm presented in Eq. (6).

Split Bregman Reconstruction Algorithm

Unconstrained convex minimization problem for CT image reconstruction models. $\text{min}{\Vert f\Vert }_{\text{TV}}+\frac{\mu }{2}{\Vert Af-p\Vert }_2^2,$ (19) where ${\Vert f\Vert }_{\text{TV}}=\left|\nabla f\right|={\displaystyle {\sum }_{s=1}^S{\displaystyle {\sum }_{v=1}^V\nabla f_{s,v}}}$

The main steps for solving the image reconstruction problem with the split Bregman algorithm for TV minimization are expressed by Eq. (19. First, the two variables are imported to transform Eq. (19) into $\begin{array}{c}\underset{f,d_x,d_y,d_z}{\min }{\left(d_x,d_y,d_z\right)}_2+\frac{u}{2}Af-p_2^2,s.t.\\ d_x={\nabla }_{x}f, d_y={\nabla }_{y}f,d_z={\nabla }_{z}f\end{array}$ (20) where ${\left(d_x,d_y,d_z\right)}_2={\displaystyle {\sum }_{s,v}\sqrt{d_{x,s,v}^2+d_{y,s,v}^2+d_{z,s,v}^2}}$. The penalty term is subsequently introduced, and the above equation is converted into $\begin{array}{c}\underset{f,d_x,d_y,d_z}{\min }{\left(d_x,d_y,d_z\right)}_2+\frac{\mu }{2}Af-p_2^2\\ +\frac{\lambda }{2}{d}_x-{\nabla }_x{f}_2^2+\\ \frac{\lambda }{2}{d}_y-{\nabla }_y{f}_2^2+\frac{\lambda }{2}{d}_z-{\nabla }_z{f}_2^2.\end{array}$ (21)

Therefore, the Bregman iteration method can be used to solve Eq. (20). $\begin{array}{l}\underset{f,d_x,d_y,d_z}{\min }{\left(d_x,d_y,d_z\right)}_2+\frac{\mu }{2}Af-p_2^2\\ +\frac{\lambda }{2}{d}_x-{\nabla }_{x}f-b_x^k{}_2^2+\\ \frac{\lambda }{2}{d}_y-{\nabla }_{y}f-b_y^k{}_2^2+\frac{\lambda }{2}{d}_z-{\nabla }_{z}f-b_z^k{}_2^2\\ b_x^{k+1}=b_x^k+\left({\nabla }_x{f}^{k+1}-d_x^{k+1}\right)\\ b_y^{k+1}=b_y^k+\left({\nabla }_y{f}^{k+1}-d_y^{k+1}\right)\\ b_z^{k+1}=b_z^k+\left({\nabla }_z{f}^{k+1}-d_z^{k+1}\right)\end{array}$ (22)

Eq. (20) can usually be updated in two alternating iterative steps f and d. $\begin{array}{ll}Step1:\hfill & f^{k+1}=\text{arg}\underset{f}{\min }\frac{u}{2}{\Vert Af-p\Vert }_2^2+\frac{\lambda }{2}{\Vert d_x^k-{\nabla }_{x}f-b_x^k\Vert }_2^2\hfill \\ \hfill & +\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\hfill \end{array}$ $\begin{array}{l}Step2:\left(d_x^{k+1},d_y^{k+1},d_z^{k+1}\right) =\text{arg}\underset{d_x,d_y,d_z}{\min }{\Vert \left(d_x,d_y,d_z\right)\Vert }_2\\ +\frac{\lambda }{2}{\Vert (d_x-{\nabla }_{x}f-b_x^k\Vert }_2^2+\frac{\lambda }{2}{\Vert (d_y-{\nabla }_{y}f-b_y^k\Vert }_2^2\\ +\frac{\lambda }{2}{\Vert (d_z-{\nabla }_{z}f-b_z^k\Vert }_2^2\end{array}$

The first step considers the specificity of the CT reconstruction matrix A and is solved using the gradient descent method. $f^{k+1}=f^k-{\alpha }_k{g}^k,$ (23) $\begin{array}{c}{g}^k=\mu ·A^T\left(A{f}^k-p\right)-\lambda ·{\nabla }_x^T\left(d_x-{\nabla }_{x}f-b_x^k\right)\\ -\lambda ·{\nabla }_y^T\left(d_y-{\nabla }_{y}f-b_y^k\right)-\lambda ·{\nabla }_z^T\left(d_z-{\nabla }_{z}f-b_z^k\right),\end{array}$ where ${\alpha }_k=\frac{{\left(g^k\right)}^{\prime }·\left(g^k\right)}{{\left(g^k\right)}^{\prime }·{\nabla }_f\left(g^k\right)·\left(g^k\right)}$ is the step size.

The second step can be solved by the generalized contraction formula. $d_x^{k+1}=\text{max}\left(s^k-\frac{1}{\lambda },0\right)\frac{{\nabla }_x{f}^k+b_x^k}{s^k},$ $d_y^{k+1}=\text{max}\left(s^k-\frac{1}{\lambda },0\right)\frac{{\nabla }_y{f}^k+b_y^k}{s^k},$ $d_z^{k+1}=\text{max}\left(s^k-\frac{1}{\lambda },0\right)\frac{{\nabla }_z{f}^k+b_z^k}{s^k},$ where $s^k=\sqrt{\left|{\nabla }_x{f}^k+b_x^k|+|{\nabla }_y{f}^k+b_y^k|+|{\nabla }_z{f}^k+b_z^k\right|}.$

The specific process of CT reconstruction based on the split Bregman algorithm is as follows.

Step 1: initial $k=0,d_x^0=d_y^0=d_z^0=b_x^0=b_y^0=b_z^0=0,f^{\left(0\right)}=0$

Step 2: Perform the gradient descent method using Eq. (21)

Step 3: Calculating intermediate d_x, d_y, d_z, b_x, b_y, b_z

Step 4: Return to step 2 until the stop condition is reached.

OS-SART

The SART is an improved iterative reconstruction algorithm proposed by Anderson and Kak in 1984. The SART is a modification of the ART in that it uses the error of all rays passing through a pixel at a certain angle to correct its value. This formula is expressed as follows: $f_j^{\left(k+1\right)}=f_j^k+\begin{array}{c}\raisebox{1ex}{${\displaystyle {\sum }_{p_i\in p_{\phi }}\left\{{\lambda }_k\frac{p_i-{\displaystyle \sum _{n=0}^{J-1}{w}_{in}{f}_{in}^{\left(k\right)}}}{{\displaystyle {\sum }_{n=0}^{J-1}{w}_{in}}}{w}_{ij}\right\}}$}\!\left/ \!\raisebox{-1ex}{${\displaystyle {\sum }_{p_i\in p_{\phi }}{w}_{ij}}$}\right.\end{array}.$

The OS-SART algorithm is an improvement over the SART algorithm that divides the projection angle using the following equation: $\begin{array}{c}{f}_j^{\left(k+1\right)}=f_j^k+\raisebox{1ex}{${{\displaystyle \sum }}_{i\in S_t}\left\{{\lambda }_k\frac{p_i-{{\displaystyle \sum }}_{n=0}^{J-1}{w}_{in}{f}_{in}^{\left(k\right)}}{{{\displaystyle \sum }}_{n=0}^{J-1}{w}_{in}}{w}_{ij}\right\}$}\!\left/ \!\raisebox{-1ex}{${{\displaystyle \sum }}_{i\in S_t}{w}_{ij}$}\right.\end{array}.$

Proposed algorithm

Based on the OS-SART algorithm and split Bregman method, we propose the 3D reconstruction algorithm OS-SART-SBTV for sparse-view neutron CT. The algorithm employs a two-step iteration process with OS-SART for image reconstruction and the split Bregman method for total variation denoising (Table 1).

Quantitative Evaluation Index

Image evaluation metrics play a vital role in assessing the quality and accuracy of reconstructed images. These metrics provide objective measures to compare different algorithms and determine their performance. The following four image evaluation metrics were selected to analyze the performance of each algorithm.

The smaller the root mean square error (RMSE), the better, and the smaller it is also implies that the closer the two images are, the more detail is retained in the original image. RMSE is defined by the following equation: $\begin{array}{c}\text{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{n=1}^N{\left({\widehat{p}}_n-p_n\right)}^2}{N}}.\end{array}$

The universal quality image (UQI) is a metric for evaluating image quality that quantifies the similarity between two images using the following equation: $\begin{array}{c}\text{UQI}=\frac{2cov\left(\widehat{\mu },\mu \right)}{{\widehat{\sigma }}^2+{\sigma }^2}·\frac{2\widehat{\mu }\mu }{{\widehat{\mu }}^2+{\mu }^2}.\end{array}$

Correlation coefficient (CC) is an evaluation metric used to quantify the correlation between two images using the following equation: $\text{CC}=\frac{{\displaystyle {\sum }_{k=1}^M\left(u\left(k\right)-\overline{u}\right)}\left(u_{\text{true}}\left(k\right)-\overline{u_{\text{true}}}\right)}{{\left\{{{\displaystyle {\sum }_{k=1}^M\left(u\left(k\right)-\overline{u}\right)}}^2{{\displaystyle {\sum }_{k=1}^M\left(u_{\text{true}}\left(k\right)-\overline{u_{\text{true}}}\right)}}^2\right\}}^{1/2}}.$

Mean structural similarity (MSSIM) is based on the concept of structural similarity and aims to measure the degree of structural similarity and distortion between two images using the following equation: $\text{MSSIM}=\frac{1}{M}{\displaystyle \sum _{i=1}^M{\left[l\left(x_i,y_i\right)\right]}^{\alpha }\times {\left[c\left(x_i,y_i\right)\right]}^{\beta }\times {\left[s\left(x_i,y_i\right)\right]}^{\gamma }},$ where brightness information is $l\left(x,y\right)=\frac{2u_x{u}_y+C_1}{u_x^2+u_y^2+C_1}$, contrast information is $c\left(x,y\right)=\frac{2{\sigma }_x{\sigma }_y+C_2}{{\sigma }_x^2+{\sigma }_y^2+C_2}$ and structure Information is $s\left(x,y\right)=\frac{{\sigma }_{xy}+C_3}{{\sigma }_x+{\sigma }_y+C_3}$.

3D Digital Head Model Experiment

In this study, the performance differences among the five algorithms were comparatively analyzed using a 3D digital head model with dimensions of 256 × 256 × 256 (Fig.1). The four iterative algorithms must be adjusted with the parameters to obtain the best-quality image. Therefore, the parameters of each algorithm were set as follows: OS-SART-TV algorithm, ${\text{λ}}_{\text{OS}-\text{SART}}=2.2$, ${\text{Niter}}_{0\text{S}-\text{SART}}=50$, ${\text{λ}}_{\text{TV}}=200$, and ${\text{Niter}}_{\text{TV}}=150$. As for the OS-SART-FGPTV algorithm, ${\text{λ}}_{\text{OS}-\text{SART}}=2.5$, ${\text{Niter}}_{\text{OS}-\text{SART}}=50$, ${\text{λ}}_{\text{FGPTV}}=0.006$, and ${\text{Niter}}_{\text{FGPTV}}=150$. As for the OS-SART-LRMA algorithm, ${\text{λ}}_{\text{OS}-\text{SART}}=2.3$, ${\text{Niter}}_{\text{OS}-\text{SART}}=50$, $\text{blcksize}=［121212］$. For the OS-SART-SBTV algorithm, ${\lambda }_{\text{OS}-\text{SART}}=2.2$, ${\text{Niter}}_{\text{OS}-\text{SART}}=50$, ${\text{λ}}_{\text{SBTV}}=0.0032$, and ${\text{Niter}}_{\text{SBTV}}=150$.

Fig. 1Full-screen View

Slice of the head model at z = 90

Figure 2 displays the dependency between the RMSE and iterations for different numbers of projection views reconstructed using the OS-SART-SBTV algorithm. The results display that OS-SART-SBTV can reach a convergence state after certain iterations. In addition, the convergence speed of the algorithm increased with an increase in sparse views, indicating that the OS-SART-SBTV algorithm can minimize the objective function and obtain a satisfactory solution under different sparse view conditions.

Fig. 2Full-screen View

Relationship between RMSE and the number of iterations of OS-SART-SBTV: (a) 30, (b) 60, and (c) 60+

Figure 3 shows the five algorithms used to reconstruct the images with different projection angles. The numbers of projection views were 10, 30, 60, and 60+. 60+ denotes the reconstruction results of the projection data after adding the Poisson and Gaussian noise models. The quality of the neutron source and the inherent noise of the system equipment can affect the acquired projection data. Therefore, by introducing Gaussian and Poisson noise to the digital head model, we can test the performance and stability of the algorithm. To analyze the noise reduction performance of the algorithm, two noise distributions were added to the digital model. The noise model parameters are the incident photon flux 1×10⁵, $m_{\text{ic}}=0$, and ${\delta }_{\text{ic}}^2=0.25$ [32]. As depicted in Figure 3, the sharpness and smoothness of the images reconstructed using the five algorithms increased progressively with respect to the number of projected views when reconstructing the noiseless projection data. When the reconstruction angle was 10 degrees, the reconstructed image quality of all five algorithms was poor, leading to an inability to differentiate the performance of the algorithms effectively. The reconstructed results of OS-SART-TV depicted in Figure 3, were superior to those of the FBP algorithm; however, its reconstructed image exhibited the problems of detail loss and excessive smoothing of edge structure, especially when reconstructing 30 views, the reconstructed image of OS-SART-TV algorithm was blurred. The third to fifth rows in Figure 3 represent the reconstructed images from the OS-SART-FGPTV, OS-SART-LRMA, and OS-SART-SBTV algorithms. The number of iterations for the four algorithms used to reconstruct the different projection views was set to 50. Specifically, the reconstruction results of the OS-SART-FGPTV and OS-SART-LRMA algorithms displayed certain improvements compared to the FBP algorithm. The OS-SART-FGPTV and OS-SART-LRMA algorithms reconstruct more details than OS-SART-TV. Compared to OS-SART-TV, OS-SART-FGPTV, and OS-SART-LRMA, the OS-SART-SBTV algorithm had fewer image artifacts and sharper edge structures when the number of reconstruction angles was 30. Compared with OS-SART-TV, the OS-SART-SBTV algorithm not only suppressed the artifacts well but also reconstructed more structural information about the image. In addition, we zoomed in on the regions of interest (ROI) for the different algorithm-reconstructed images depicted in Fig. 3. According to the ROI of the FBP algorithm, serious artifacts were present in the ROI of the reconstructed image, and a significant amount of detailed image information was lost. It is well known that the ROI edges of images reconstructed by OS-SART-TV are excessively smooth, resulting in blurred images. Compared with the OS-SART-FGPTV and OS-SART-LRMA algorithms, the OS-SART-SBTV algorithm displayed fewer artifacts in the ROI images and reconstructed more detailed image information and edge structure features.

Fig. 3Full-screen View

Slice of reconstructed images of the head model obtained using FBP, OS-SART-TV, OS-SART-FGPTV, OS-SART-LRMA, and OS-SART-SBTV at z = 90

Figure 4 represents the relationship between the number of iterations and the variation in the RMSE for the different algorithms. According to the curves in Fig. 4, the OS-SART-SBTV algorithm has the fastest convergence rate, and the RMSE values for each iteration of the OS-SART-SBTV algorithm were smaller than those of the other two algorithms.

Fig. 4Full-screen View

Relationship between RMSE and the iterations for different algorithms: (a) 30, (b) 60, and (c) 60+

Similar conclusions were drawn from Fig. 5, which displays the error images for each reconstruction algorithm. As shown in Fig. 5, different algorithms reconstructed the error images for the 30, 60, and 60+ projection views. Although we could not discriminate the superiority of each algorithm when there were 60 views, the OS-SART-SBTV algorithm showed great superiority when the number of projected views was small, for example, 30 views. From the error image in Fig. 5, the error image of the FBP had serious artifacts and lost a lot of image structure information, indicating that the FBP algorithm could not be applied to sparse reconstruction. Compared to FBP, the four iterative algorithms performed well in reconstructing detailed image information and suppressing artifacts. The OS-SART-SBTV algorithm exhibited the least loss of detailed structural information in the error image. To further compare and analyze the performance differences between the algorithms, we scaled up the ROI of the error images of the different algorithms (red rectangular area in Fig. 5). According to the ROI of the error image of FBP, it is evident that several image structure features were lost, and multiple artifacts appeared in the error image. By comparing the ROI of the OS-SART-TV, OS-SART-FGPTV, OS-SART-LRMA, and OS-SART-SBTV error images, the OS-SART-SBTV algorithm could reconstruct more detailed structures. Therefore, the OS-SART-SBTV algorithm outperformed the other two iterative algorithms in terms of reconstructing the image detail structures.

Fig. 5Full-screen View

Slice of error images at z = 90

Figure 6 displays the vertical and horizontal profiles of the reconstructed FBP, OS-SART-TV, OS-SART-FGPTV, OS-SART-LRMA, and OS-SART-SBTV. The results resemble the reconstructed and reference images. As shown in Fig. 6, fluctuations were present in the profile lines of FBP, and the pixel values of the reconstructed image significantly deviated from the actual values. Although the four iterative algorithms were superior to FBP, they deviated from the real pixel values. The OS-SART-SBTV algorithm reconstructs images with pixel values closest to the real pixel values. Compared to the FBP, OS-SART-TV, OS-ART-FGPTV, and OS-SART-LRMA algorithms, the reconstructed images of the OS-SART-SBTV algorithm were closer to the reference image, regardless of the vertical or horizontal profiles. This illustrates the excellent performance of the OS-SART-SBTV algorithm for sparse-view CT 3D reconstruction.

Fig. 6Full-screen View

Profile line graphs of reconstructed images by different algorithms: vertical and horizontal

Based on the above results, the FBP algorithm cannot be applied to sparse reconstructions. For further comparison, we selected the RMSE, CC, MSSIM, and UQI to illustrate the edge protection and noise suppression effects of the four iterative algorithms.

Figure 7 depicts the UQI, RMSE, CC, and MSSIM evaluation metrics for the reconstructed images using the four algorithms. First, we quantitatively analyzed images reconstructed using noiseless projection data. Figure 7a displays that the RMSE value of the algorithm increased as the number of projections decreased. When the algorithms reconstruct the projection angle in the same manner, the RMSE of the OS-SART-SBTV reconstructed image was minimized. When the reconstructed sparse view was 30, the RMSE of OS-SART-SBTV was 0.0246. As illustrated in Fig. 7, the CC, MSSIM, and UQI of all algorithms decreased with the number of projection angles. When the number of reconstructed angles was the same, the CC, MSSIM, and UQI of OS-SART-SBTV were at the maximum, indicating that the algorithm reconstructed the image with the best quality. We subsequently quantitatively analyzed the reconstructed images using the noise-containing projection data. The RMSE analysis revealed that the OS-SART-SBTV reconstructed image had the lowest RMSE value when reconstructing the same amount of noise-containing projection data. From the CC, MSSIM, and UQI, the OS-SART-SBTV reconstructed images had the greatest values when reconstructing the same number of noise-containing projection data. In conclusion, according to the reconstruction results of each algorithm, the image quality of our proposed algorithm was the highest when reconstructing the same number of projection views.

Fig. 7Full-screen View

Image evaluation metrics of different algorithms for reconstructing head model. (a) RMSE, (b) CC, (c) MSSIM, and (d) UQI

Real Neutron Beam Projection Experiment

We analyzed the performance of the OS-SART-SBTV algorithm in real applications by reconstructing the projection data of the clock model (Fig.8). The clock model-based neutron CT projection data were supplied by Schillinger [33]. Schillinger presented 201 equirectangular neutron photographs over a 180-degree range.

Fig. 8Full-screen View

Slice of clock model at z = 110

We used five algorithms to reconstruct the clock-projection data from a real neutron CT experiment to verify the superiority of our algorithms. As shown in Fig. 9, the reconstructed image of FBP contains obvious artifacts. The reconstructed image of the OS-SART-TV was extremely smooth, causing a relatively serious loss of detailed information. Although the OS-SART-FGPTV and OS-SART-LRMA algorithms improved the reconstructed images compared to the OS-SART-TV algorithms, certain structural information was still lost in the images reconstructed by the OS-SART-FGPTV and OS-SART-LRMA algorithms when the sparsity was very small. Visual analysis revealed that the OS-SART-SBTV was better than the other algorithms at eliminating artifacts and suppressing noise. In addition, the OS-SART-SBTV algorithm preserved more detailed structural information, and the reconstructed images were of the highest quality.

Fig. 9Full-screen View

Slices of reconstructed images at z = 110 for different algorithms

The images of the real neutron projection data reconstructed by the four algorithms were evaluated using four image evaluation metrics (Fig. 10). The RMSE of OS-SART-SBTV was the lowest, indicating that the algorithm reconstructed the image nearest to the original image. OS-SART-SBTV had the maximum UQI, CC, and MSSIM, indicating the best performance of the algorithm. In summary, according to the reconstruction results, OS-SART-SBTV is the best method for removing noise, suppressing artifacts, and reconstructing detailed structures.

Fig. 10Full-screen View

Image evaluation metrics of different algorithms for reconstructing clock model. (a) RMSE, (b) CC, (c) MSSIM, and (d) UQI