Multiobjective optimization guided by image quality index for limited-angle CT image reconstruction

Abstract

BACKGROUND:

Due to the incomplete projection data collected by limited-angle computed tomography (CT), severe artifacts are present in the reconstructed image. Classical regularization methods such as total variation (TV) minimization, ℓ₀ minimization, are unable to suppress artifacts at the edges perfectly. Most existing regularization methods are single-objective optimization approaches, stemming from scalarization methods for multiobjective optimization problems (MOP).

OBJECTIVE:

To further suppress the artifacts and effectively preserve the edge structures of the reconstructed image.

METHOD:

This study presents a multiobjective optimization model incorporates both data fidelity term and ℓ₀-norm of the image gradient as objective functions. It employs an iterative approach different from traditional scalarization methods, using the maximization of structural similarity (SSIM) values to guide optimization rather than minimizing the objective function.The iterative method involves two steps, firstly, simultaneous algebraic reconstruction technique (SART) optimizes the data fidelity term using SSIM and the Simulated Annealing (SA) algorithm for guidance. The degradation solution is accepted in the form of probability, and guided image filtering (GIF) is introduced to further preserve the image edge when the degradation solution is rejected. Secondly, the result from the first step is integrated into the second objective function as a constraint, we use ℓ₀ minimization to optimize ℓ₀-norm of the image gradient, and the SSIM, SA algorithm and GIF are introduced to guide optimization process by improving SSIM value like the first step.

RESULTS:

With visual inspection, the peak signal-to-noise ratio (PSNR), root mean square error (RMSE), and SSIM values indicate that our approach outperforms other traditional methods.

CONCLUSIONS:

The experiments demonstrate the effectiveness of our method and its superiority over other classical methods in artifact suppression and edge detail restoration.

Keywords

CT reconstruction limited-angle CT image quality assessment multiobjective optimization

1 Introduction

Limited-angle X-ray CT is restricted by the scanning object or scanning environment, such as the imaging of breast CT [1], C-arm CT [2] and dental CT [3]. Typically, it can only scan the object within a limited scanning range, resulting in incomplete collection projection data.In limited-angle CT image reconstruction problems, the images reconstructed by filtered back projection (FBP) [4, 5], which will have obvious artifacts near the edges. For some traditional iterative image reconstruction algorithms [6], such as Simultaneous Algebraic Reconstruction (SART) [7] and Expectation Maxmization (EM) [8], it is challenging to reconstruct images with clear edges under conditions of incomplete projection data.

To solve the problems of artifacts and detail loss in reconstructed images. Regularization strategie has been proposed for image reconstruction of limited-angle CT. A regularization strategie is to reduce the ill-posed problem by adding prior information. TV minimization [9, 10] use the ℓ₁-norm of the amplitude of image gradient as the regularization term, can effectively preserve the edge of the image. However, for limited-angle reconstruction problem, the reconstruction image may be oversmooth and there will still be some distortion near the edge with TV minimization. To further suppress artifacts and preserve the edges, some variants of total variation [11 –16] have demonstrated that they can to mitigate the over-smoothness of reconstructed images to a certain extent. For example, an adaptive-weighted TV (AwTV) [14] algorithm is proposed, which considers the anisotropic edge property. The AwTV algorithm can adaptively adjust the weights according to the local image-intensity gradient and the edge details can be effectively preserved. The iteratively reweighted anisotropic TV (ATV) method [15] which combines the reweighted technique with ATV, can effectively utilize the property of image sparsity. The local anisotropic total variation repair method of sinogram (LATV) [16] leverages the continuity of local gradients in the sinogram to improve local high-frequency characteristics that are blurred by conventional TV algorithm. However, although these variants mitigate the problem of distortion, the reconstruction image will still be oversmooth. To emphasize the sparsity of image under the gradient transform and preserve the edge of the image while suppressing artifacts, the ℓ₀-norm of the image gradient is taken as a regularization term [17], the ℓ₀ minimization of image gradient was proposed, which can effectively preserve the image edge, but some subtle details of image will still be lost. To further suppress artifacts and restore image edge details, some regularization reconstruction methods utilize wavelet frame to suppress the artifacts [18 –22], as models based on wavelet frame can sparsely approximate the piecewise smooth function like images. For example, combining wavelet frame with ℓ₀ minimization, where the ℓ₀-norm of wavelet coefficients is penalized. Guided image filtering (GIF) [23] is a smoothing operator with edge preserving, which is effective in many applications,including edge-aware smoothing, detail enhancement, joint upsampling. To further improving the reconstructed image quality, the GIF was introduced into the reconstruction process [24 –26]. The method that couples ℓ₀ minimization with GIF [25] can effectively preserve the image edge while smoothing the image. And by leveraging the advantages of GIF, the method [26] transfers reconstructed results from the TV method to the ℓ₀-norm of a gradient image-based total variation minimization (TVL0) method. This method better preserves edge structures and suppresses limited-angle artifacts and noise.

To solve most regularization reconstruction models as above are usually utilize the following two optimization problems [9 , 27–29]: $min_{f \geq 0} {\frac{1}{2} {∥ Af - g^{δ} ∥}^{2} + aR (f)},$ (1) or $min_{f \geq 0} R (f) s . t . \frac{1}{2} {∥ Af - g^{δ} ∥}^{2} \leq ɛ,$ (2) where ∥Af - g^δ ∥ ² represents a data fidelity term, which is used to control the data errors, a > 0 is a regularization parameter, R (f) represents a regularization term, and ɛ is an upper bound for the objective function. The problems Equations (1) and (2) are the weighting method and ɛ-constraint method in the scalarization methods [30, 31] of the following multiobjective optimization problem (MOP) respectively: $min_{f \geq 0} {\frac{1}{2} {∥ Af - g^{δ} ∥}^{2}, R (f)} .$ (3) However, there are many weaknesses in solving MOP Equation (3) with scalarization methods. Using the weighting method requires assigning different weights to different objective functions. Although it is simple and convenient, it is difficult to set appropriate weights to obtain a Pareto optimal solution. Moreover, if the Pareto-optimal front of the MOP is non-convex, it is difficult for us to find all Pareto optimal solutions. For the ɛ-constraint method, we need to choose the most important sub-objective function as the main objective from all objective functions, and take the remaining sub-objective functions as constraints, and set an upper bound for each sub-objective function. Although the optimal solutions of the ɛ-constraint method are Pareto optimal solutions, we need to take the optimized value of each sub-objective function as a constraint, so it is time consuming.

The scalarization method refers to a parameter-dependent single-objective optimization problem, which is proposed for MOP. For the scalarization method, the selection of the parameter is very important, and using scalarizations may destroy the structure of the optimization problem. Non-scalarization methods may have advantages for optimization problems, such as lexicographic optimality [32]. In lexicographic optimization, we need to consider the lexicographic order when comparing objective functions. We first arrange the objective functions according to their absolute importance, minimizing the objective function f_k in the sorted order and using the optimal objective values of f_i, i < k as constraints, until the single-objective optimization problem has a unique optimal solution.

In recent years, many deep learning (DL) methods have been proposed for CT image denoising and have achieved many encouraging results [33, 34]. The residual encoder-decoder convolutional neural network (RED-CNN) [35] is a combination of autoencoder, deconvolution network and shortcut connection, which has shown potential in suppressing artifacts and preserving image structure. A method based on a generative adversarial network (GAN) with Wasserstein distance and perceptual similarity has been proposed that can effectively suppress artifacts [36]. However, DL needs a lot of data to train the network, and in many applications, label data can not be obtained in large quantities and conveniently. Using a CNN or DL to solve a CT related inverse problem is no better than the traditional regularization methods [37]. Therefore, there are still many challenges in applying neural networks to CT reconstruction, so we only consider traditional optimization algorithms.

In this paper, we proposed a new method for limited-angle CT image reconstruction using a multiobjective optimization model guided by the image quality index. Unlike traditional scalarization methods that solely minimize the objective function, our method introduce the SSIM to guide the optimization process. The multiobjective optimization model considers the data fidelity term and ℓ₀-norm of the image gradient as two objective functions. Using traditional scalarization algorithms to solve the MOP aims to minimize the objective function, however, the minimum objective function value does not represent the highest image quality. For example, if the gradient is zero, the image is black. Therefore, we employ an iterative method instead of using a scalarization algorithm to solve the model. Our method is structured into two steps, and repeat these two steps in an iteration format, and the result of the previous step is used as the initial value of the next step. The first step is to optimize the data fidelity term by SART, then calculate the SSIM value between the reference image and the iterative points that are solved by SART. The optimization process will be stopped when the image quality degrades. Then, the SA algorithm [38] is introduced to accept the degradation solution in a probabilistic form to avoid getting into the local minimum. If we accept the degradation solution, continue to repeat SART, otherwise, we use the GIF method to produce an alternative solution because GIF can effectively preserve the image edges before proceeding to the second step. The second step is to take the result of the first step as a constraint, then use ℓ₀ minimization to optimize ℓ₀-norm of the image gradient, the SSIM, SA algorithm and GIF are still introduced to guide the optimization process as in thefirst step.

The remainder of this paper is organized as follows. In Section 2, the proposed multiobjective optimization model and solving algorithm are presented, along with the related methods required for our algorithm. In Section 3, simulated experiment and real data experiments are conducted to evaluate the effectiveness of the proposed method. In Section 4, the conclusions are presented.

2 Model and method

The discrete reconstruction problem of limited-angle CT is formulated as: $g^{δ} = Af,$ (4) where f ∈ R^N is the image to be reconstructed, g^δ ∈ R^M is the projection data containing noise, A = (a_ij) _M×N ∈ R^M×N (M ⪡ N) is the system matrix, a_ij is the weight and i is the X-ray index, and j is the pixel index of the reconstructed image. Due to the incompleteness of the projection data collected, the reconstruction result of the classical algebraic reconstruction algorithm is unsatisfactory, such as SART. Therefore, we proposed our method to solve the limited-angle CT reconstruction problem.

2.1 Guided image filtering

The guided image filtering [23] is an explicit image filter, which can be used to preserve image edges, suppress artifacts and enhance details. The method is shown as follows: $f^{out} = GIF (I, f^{input}) .$ (5) where GIF is a function of guided image filtering, f^input is the filtering input image, I is the guidance image and f^out is the filtering output image.

The key to defining guide filter lies in assuming a local linear model between the guidance image I and the filtering output image. Within the square window ω_k of radius r of the pixel k, the linear transform can be formulated as: $f_{i}^{out} = a_{k} I_{i} + b_{k}, \forall i \in ω_{k},$ (6) where (a_k, b_k) are some constant linear coefficients in ω_k. The difference between the input image f_input and the filtering output image f_out can be expressed as the input image minus some unnecessary components like noise. By minimizing the following the cost function, an optimal solution of a_k and b_k can be obtained: $E (a_{k}, b_{k}) = \sum_{i \in ω_{k}} ({(a_{k} I_{i} + b_{k} - f_{i}^{input})}^{2} + ɛ a_{k}^{2}),$ (7) where ɛ is a regularization parameter used to penalize large a_k. The solution can be expressed as: $a_{k} = \frac{\frac{1}{| ω |} \sum_{i \in ω_{k}} I_{i} f_{i}^{input} - μ_{k} {\bar{f}}_{k}^{input}}{σ_{k}^{2}}, b_{k} = {\bar{f}}_{k}^{input} - a_{k} μ_{k},$ (8) where |ω| is the number of pixels in ω_k, μ_k is the mean of I within ω_k, $σ_{k}^{2}$ is the variance of I within ω_k, and ${\bar{f}}_{k}^{input} = \frac{1}{| ω |} \sum_{i \in ω_{k}} f_{i}^{input}$ denotes the mean of $f_{i}^{input}$ in ω_k. Due to one pixel may belong to different windows, resulting in different filtering outputs, the final output is obtained by averaging across all overlapping windows containing the pixel. Thus, the filtering output can formulate as: $f_{i}^{out} = \frac{1}{| ω |} \sum_{k | i \in ω_{k}} (a_{k} I_{i} + b_{k}) = {\bar{a}}_{i} I_{i} + {\bar{b}}_{i} .$ (9) Here, denoting ∑_{k|i∈ω_k}a_k = ∑_{k∈ω_i}a_k. Where, ${\bar{a}}_{i} = \frac{1}{| ω |} \sum_{k \in ω_{i}} a_{k}$ , ${\bar{b}}_{i} = \frac{1}{| ω |} \sum_{k \in ω_{i}} b_{k}$ are the average coefficients of all windows overlapping pixel i.

2.2 Proposed model and method

The ℓ₀-norm of the image gradient minimization can effctively suppress streak artifacts and preserve the image edges. Thus, we considered the ℓ₀-norm of the image gradient as one of the optimization objective functions. To control the data errors, we incorporate the data fidelity term as another objective function.

To further improve the quality of the reconstructed image, we proposed a multiobjective optimization model for limited-angle CT image reconstruction, which is formulated as follows: $min_{f} {\frac{1}{2} {∥ Af - g^{δ} ∥}_{D}^{2}, {∥ \nabla f ∥}_{0}}, s . t . f \geq 0,$ (10)

where ${∥ Af - g^{δ} ∥}_{D}^{2}$ is the data fidelity term, D is the diagonal matrix with diagonal elements $\frac{1}{\sum_{j = 1}^{n} a_{ij}} (\sum_{j = 1}^{n} a_{ij} \neq 0)$ , ∥∇ f ∥ ₀ = ∑_i,j # { (i, j) ||∂_xf_i,j| + |∂_yf_i,j| ≠ 0 }, #{ } is the counting operator that counts the number of pixels satisfying |∂_xf_i,j| + |∂_yf_i,j| ≠ 0. The image gradient at (i, j) is given as ∇f_i,j = (∂_xf_i,j, ∂_yf_i,j) = (f_i,j - f_i-1,j, f_i,j - f_i,j-1).

The scalarization methods are usually used to solve multiobjective optimization problems, such as weighting method and ɛ-constraint method. However, there are some shortcomings of scalarization method. For example, using weighting method to solve the MOP Equation (10) cannot find all Pareto optimal solutions. The other nonscalarizing methods, such as lexicographic optimization, require considering the MOP in a hierarchical manner. Moreover, both methods are based on minimizing objective function under constraints, but the minimum objective function does not mean the highest image quality. For example, minimum-∥ ∇ f ∥ ₀ can result in the reconstructed image being oversmooth.

During the reconstruction process, we introduce the SSIM as a key evaluation criterion for the reconstructed images. SSIM evaluates image quality by focusing on the preservation of structural information. Reconstructed images often suffer from artifacts or oversmoothing, which can obscure fine structural details. Compared to other image quality indices like Peak Signal-to-Noise Ratio (PSNR) and Root Mean Square Error (RMSE), SSIM is more sensitive to structural details. In practical applications, the clarity and structural integrity of the reconstructed image are crucial. By prioritizing structural aspects, SSIM ensures that images with higher structural quality and clearer details are produced during the optimization process. The GIF we introduced after the degradation point is rejected, one of the significant advantages of GIF is its ability to smooth images while effectively preserving edges and fine details. Traditional smoothing methods, such as mean filtering and Gaussian filtering, often blur edges and lose critical details. GIF, by using a guided image, can reduce noise while maintaining the structural integrity of important features, which is crucial for high-quality image reconstruction.

In this paper, we address solving MOP Equation (10) directly through an iterative method. During the optimization process, SSIM is introduced to guide the optimization of the objective function. When the SSIM values degrade, the SA algorithm is introduced to probabilistically accept the degraded solution to avoid getting trapped in a local minimum. Subsequently, the GIF is introduced to further improve image quality. The detailed process is as follows:

Step 1: Optimizing the data fidelity term establishes a solid foundation for further refinements. By minimizing this term first, we can ensure accuracy of data reconstruction. First, to optimize the following data fidelity term: $min_{f} \frac{1}{2} {∥ Af - g^{δ} ∥}_{D}^{2}, s . t . f \geq 0,$ (11) we solved the optimization problem Equation (11) by the following SART to obtain the initial point: $f^{(k + 1)} = f^{(k)} - ω {DA}^{T} W^{- 1} ({Af}^{(k)} - g^{δ}),$ (12) where k represents the number of iterations, ω ∈ (0, 2) is a relax factor, which controls the step size of each iteration. D is a diagonal matrix with diagonal elements $\frac{1}{\sum_{j = 1}^{n} a_{ij}} (\sum_{j = 1}^{n} a_{ij} \neq 0)$ , W is also a diagonal matrix whose diagonal elements are the sum of the weights of each projection, with the diagonal elements being $\sum_{l = 1}^{N} a_{il}$ . A and A^T represents the forward projection and the back projection respectively.

Then, we introduce the image quality index SSIM, SA algorithm and GIF to guide the optimization process. We calculate the SSIM values between reference image I and iteration point f^(k+1), f^(k) respectively (note it as SSIM (f^(k+1), I) and SSIM (f^(k), I)). If the image quality improves during the iteration process, i.e., if SSIM (f^(k+1), I) is lager than SSIM (f^(k), I), we repeat the iteration process Equation (12). Otherwise, we accept the iteration point f^(k+1) in the form of probability to avoid getting into local minimum, the form of probability is as follows: $P (f^{(k + 1)}) = {\begin{matrix} 1, if Δ E_{1} < 0 \\ e^{- \frac{Δ E_{1}}{T_{1}}}, if Δ E_{1} > 0 \end{matrix},$ (13) where ΔE₁ = SSIM (f^(k), I) - SSIM ((f^(k+1), I), T₁ is the initial temperature value of the temperature set by SA algorithm.

If accept the iteration point f^(k+1), we repeat the iteration process Equation (12), otherwise we introduce guided image filtering [23] to improve image quality, which is given as follows: ${\hat{f}}^{(k + 1)} = GIF (f^{(k)}, f^{(k + 1)}),$ (14) where GIF is a function of guided image filtering, f^(k+1) is a filtering input image, f^(k) is a guidance image, and the filtering output image is denoted as ${\hat{f}}^{(k + 1)}$ . Then, go to step 2 with ${\hat{f}}^{(k + 1)}$ as the initial point.

Step 2: Add ${\hat{f}}^{(k + 1)}$ as a constraint, which can guarantee the data fidelity function keep its optimal value, and then optimize the ℓ₀-norm of the image gradient with given point ${\hat{f}}^{(k + 1)}$ , which is given as follows: $min_{f} {∥ \nabla f ∥}_{0}, s . t . \frac{1}{2} {∥ f - {\hat{f}}^{(k + 1)} ∥}^{2} < ɛ, f \geq 0,$ (15) where ${∥ f - \hat{f^{(k + 1)}} ∥}^{2}$ is the constraint term, ∥ ∇ f ∥ ₀ is the ℓ₀-norm of the image gradient. By introducing parameter η, which combines the constraint term and objective function [10], the Equation (15) is convert into the following Equation (16): $min_{f} η {∥ \nabla f ∥}_{0} + \frac{1}{2} {∥ f - {\hat{f}}^{(k + 1)} ∥}^{2} s . t . f \geq 0,$ (16) where η > 0 is a regularization parameter. To realize the minimization of ℓ₀-norm of the image gradient, the Equation (16) is solved by ℓ₀ minimization method [17]. We introduce two auxiliary variables h and v, which are close to the corresponding gradient components in x and y direction. The problem Equation (16) can be approximately treated as the following optimization problem: $min_{f, h, v} {\sum_{i, j} {(f - {\hat{f}}^{(k + 1)})}^{2} + η^{*} C (h, v) + μ [{(\partial_{x} f_{i, j} - h_{i, j})}^{2} + {(\partial_{y} f_{i, j} - v_{i, j})}^{2}]},$ (17) where η^* is positive parameter, C (h, v) =∑_i,j { (i, j) ||h_i,j| + |v_i,j| ≠ 0 }, μ is an automatically parameter in interations to guarantees the variables (h_i,j, v_i,j) to close (∂_xf_i,j, ∂_yf_i,j).

To solve the problem Equation (17), let t be the iteration index. The problem Equation (17) is decomposed into the following two sub-problems, solved for (h, v) and f by alternating iteration method respectively: $(h^{(t + 1)}, v^{(t + 1)}) = \underset{h, v}{arg min} {\sum_{i, j} [{(\partial_{x} f_{i, j} - h_{i, j})}^{2} + {(\partial_{y} f_{i, j} - v_{i, j})}^{2}] + η^{*} C (h, v)},$ (18) $f^{(k + 2, t + 1)} = \underset{f}{arg min} {\sum_{i, j} {(f - {\hat{f}}^{(k + 1)})}^{2} + μ^{(t)} [{(\partial_{x} f_{i, j} - h_{i, j}^{(t)})}^{2} + {(\partial_{y} f_{i, j} - v_{i, j}^{(t)})}^{2}]} .$ (19)

To solve the sub-problem Equation (18), the solution can obtain by $(h_{i, j}^{(t + 1)}, v_{i, j}^{(t + 1)}) = {\tilde{H}}_{\frac{η}{μ^{(t)}}} (| h_{i, j} |, | v_{i, j} |) = {\begin{matrix} (0, 0), {(\partial_{x} f_{i, j}^{(k + 2, t + 1)})}^{2} + {(\partial_{y} f_{i, j}^{(k + 2, t + 1)})}^{2} < \frac{η}{μ^{(t)}} \\ (\partial_{x} f_{i, j}^{(k + 2, t + 1)}, \partial_{y} f_{i, j}^{(k + 2, t + 1)}), otherwise \end{matrix} .$ (20)

To solve the sub-problem Equation (19) with h^(t) and v^(t), the analytical solution can obtain by: $f^{(k + 2, t + 1)} = F^{- 1} {\frac{F ({\hat{f}}^{(k + 1)}) + μ^{(t)} (F^{*} (\partial x) F (h^{(t)}) + F^{*} (\partial y) F (v^{(t)}))}{F (1) + μ^{(t)} (F^{*} (\partial x) F (\partial x) + F^{*} (\partial y) F (\partial y))}},$ (21) where F denotes Fast Fourier transform (FFT), F^-1 denotes inverse Fourier transform, F^* () denotes complex conjugation, F (1) denotes the Fourier transform of the delta function. The iteration format of the solution of Equation (17) can be rewritten as follows: ${\begin{matrix} f^{(k + 2, t + 1)} = F^{- 1} {\frac{F ({\hat{f}}^{(k + 1)}) + μ^{(t)} (F^{*} (\partial x) F (h^{(t)}) + F^{*} (\partial y) F (v^{(t)}))}{F (1) + μ^{(t)} (F^{*} (\partial x) F (\partial x) + F^{*} (\partial y) F (\partial y))}} \\ (h_{i, j}^{(t + 1)}, v_{i, j}^{(t + 1)}) = {\tilde{H}}_{\frac{η}{μ^{(t)}}} (| h_{i, j} |, | v_{i, j} |) \\ μ^{(t + 1)} = κ μ^{(t)}, (κ > 1) \end{matrix},$ (22) where parameter μ is multipled by κ in each iteration to speed up the convergence in the process of alternative iterative solution. We can obtain f^(k+2) = f^(k+2,t_max), where the t_max is the maximum iteration number. Then, we also introduce image quality index SSIM, SA algorithm and GIF to avoid the image quality degradation. If SSIM (f^(k+2), I) lager than $SSIM ({\hat{f}}^{(k + 1)}, I)$ , we repeat the iteration process Equation (22), otherwise SA algorithm is introduced to avoid falling into local optimum and accept iteration point f^(k+2) in the form of probability as follows: $P (f^{(k + 2)}) = {\begin{matrix} 1, if Δ E_{2} < 0 \\ e^{- \frac{Δ E_{2}}{T_{2}}}, if Δ E_{2} > 0 \end{matrix},$ (23) where $Δ E_{2} = SSIM ({\hat{f}}^{(k + 1)}, I) - SSIM (f^{(k + 2)}, I)$ and T₂ is the initial temperature value of SA algorithm. If the point of f^(k+2) is accepted, we repeat the iteration process Equation (22). Otherwise we use guided image filtering to get a new alternative point ${\hat{f}}^{(k + 2)} = GIF ({\hat{f}}^{(k + 1)}, f^{(k + 2)}),$ (24) where f^(k+2) is a filtering input image, ${\hat{f}}^{(k + 1)}$ is a guidance image, ${\hat{f}}^{(k + 2)}$ is a filtering output image. Then, go to step 1 with ${\hat{f}}^{(k + 2)}$ , which is used as the initial point.

We repeat these two steps in iteration format, and take the result of the previous step as the initial point of the next step in the iterative process until the stop conditions are met. The proposed method minimize one of the multiobjective function according to maximize the image quality index (abbr. MM), and the SA algorithm is introduced to avoid falling into a local minimum by accepting the degradation solution probabilistically. If the degradation solution is rejected, stop the optimization process and use the GIF method to produce an alternative solution, which can further improve the image quality, as GIF can effectively protect image edges. The proposed method is named as MM-SA-GIF method and will be shown in Algorithm 1.

Algorithm 1 The proposed MM-SA-GIF method

3 Simulated and real data experiments

In this section, simulated data and real projection data are utilized to evaluate the performance of the proposed MM-SA-GIF method in suppressing artifacts and preserving image edges. To validate the effectiveness of the proposed method, the classical iterative reconstruction methods TV minimization [10], ℓ₀ minimization [17], ℓ₀-GIF method [25] and lexicographic optimality of TV method (LO-TV) are used for comparison. All parameters are chosen to obtain a higher quality reconstruction image by trial and error.

3.1 Simulated NCAT phantom experimental results

To evaluate and validate our proposed algorithm, a NCAT (NURBS-based Cardiac-Torso) phantom with size of 512 is used in simulated experiment, Fig. 1 shows the NCAT phantom and the photon count for poisson noise is 1.0 × 10⁵. The geometry scanning parameters for simulated CT imaging system are listed in Table 1.

Fig. 1

NCAT phantom.

Table 1

The geometrical scanning parameters

Parameter	Value
Distance from source to rotation center	800mm
Distance from detector to rotation center	400mm
Total number of detectors	1024
The interval of projection view	1°
Pixel size of reconstruction images	1mm
Size of reconstructed image	512×512

Three scanning ranges [0, 120°] and [0, 150°] are investigated. The reconstruction parameters for each algorithm are shwon as follows: for TV minimization: the parametes, N_TV = 20, α_TV = 0.2; For ℓ₀ minimization: (1) for the scanning range [0, 120°], the smoothing parameter η = 0.0025, κ = 15; (2) for the scanning range [0, 150°], η = 0.003, κ = 15. Among the above two methods, the stopping condition is that reaches the number of iterations N = 5000 and 4000 for the scanning ranges [0, 120°] and [0, 150°] respectively. For ℓ₀ - GIF method: (1) for the scanning range [0, 120°], N = 3500, η = 0.007, κ = 15, r = 4, ɛ = 1.1 × 10^-5; (2) for the scanning range [0, 150°], N = 2500, η = 0.006, κ = 15, r = 4, ɛ = 1.1 × 10^-5. For LO-TV method: (1) for the scanning range [0, 120°], N_SART = 500, N_TV = 300, α_TV = 0.0045; (2) for the scanning range [0, 150°], N_SART = 500, N_TV = 300, α_TV = 0.005; For our method MM-SA-GIF, the reconstruction parameters are given in the following Table 2. It is important to note that for our proposed method, the set number of iterations is the maximum number of iterations. During the optimization process, an interruption mechanism will be employed to terminate the optimization early when the image quality no longer improves.

Table 2

The reconstruction parameters of MM-SA-GIF of NCAT Experiment

Scanning range	Parameters
[0, 120°]	J = 400, K₁ = 10, K₂ = 200, η = 0.005, κ = 15, r = 2, ɛ = 1.1 × 10^-5, T₁ = 20, T₂ = 20
[0, 150°]	J = 300, K₁ = 10, K₂ = 200, η = 0.007, κ = 15, r = 2, ɛ = 1.1 × 10^-5, T₁ = 100, T₂ = 100

Figures 2 and 3 show the reconstruction images by TV minimization, ℓ₀ minimization, ℓ₀-GIF method and our proposed method for the scanning ranges [0, 120°] and [0, 150°]. From Fig. 2, it can be observed that the reconstruction image by TV minimization cannot fully restore the edges and fails to suppress artifacts.The image reconstructed by ℓ₀ minimization exhibits a more complete structure compared to TV minimization, but it still fails to suppress artifacts, as seen in the rectangle in the lower left corner. Although the image reconstructed by the ℓ₀-GIF method suppresses the artifacts in the lower left corner, it does not suppress the artifacts indicated by the arrow in the middle. The proposed MM-SA-GIF method reconstructs images with significantly more complete and clearer edges compared to other methods and effectively suppresses artifacts. From Fig. 3, it can be observed that as the projection data increases, the artifacts in the images reconstructed by TV minimization decrease, but some artifacts still remain. Although ℓ₀ minimization suppresses most of the artifacts, the edge structures are still not clear enough. Compared to ℓ₀ minimization, the ℓ₀-GIF method struggles to balance artifact suppression and over-smoothing during the experiment. As a result, while the ℓ₀-GIF images appear smoother, some artifacts still remain. From the enlarged detail images, it can be seen that the MM-SA-GIF method presents clearer edge structures and preserves contour integrity.

Fig. 2

The results of simulated NCAT experiment reconstructed images using different methods for the scanning range [0, 120°].The top left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the bottom left corners of each image represent enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored. The display window for reconstruction image is [0,1].

Fig. 3

The results of simulated NCAT experiment reconstructed images using different methods for the scanning range [0, 150°]. The top left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the top left and bottom right corners of each image represent enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed images are restored. The display window for reconstruction image is [0,1].

Figure 4 shows the horizontal profiles (300th row) of the reconstructed image of simulated NCAT experiment by different algorithms from different scanning ranges of [0, 120°], [0, 150°], which can evluate the quality of the reconstructed image more intuitively. As can be seen from Fig. 4, the horizontal profiles obtained by MM-SA-GIF method are colser to the reference image, indicating that our proposed method retores more detailed information of the image.

Fig. 4

Horizontal profiles (300th row) of the reconstructed images of simulated NCAT experiment by different algorithms in different scanning ranges [0, 120°], [0, 150°].

To assess the quality of reconstructed images more directly, we introduced image quality metrics such as peak signal-to-noise ratio (PSNR), Root Mean Square Error (RMSE), and SSIM to evaluate the performance of different reconstruction methods. From Table 3, it can be observed that when the scanning ranges are [0, 120°] and [0, 150°], the PSNR and RMSE values of images reconstructed by the MM-SA-GIF method are significantly superior to the other three methods. We need to acknowledge that the SSIM values of the other three methods may be higher than those of our algorithm. This is because our algorithm utilizes the SA algorithm, which probabilistically accepts degraded solutions to explore a larger solution space, resulting in lower SSIM values compared to other algorithms. However, as evident from the reconstructed images, our algorithm significantly outperforms others in suppressing artifacts and restoring structures.

Table 3

The image quality index values of NCAT phantom experiment reconstruction iamge

Scanning range	Algorithm	PSNR	RMSE	SSIM
[0, 120°]	TV	30.4475	0.0310	0.9735
	ℓ₀	31.0019	0.0291	0.9690
	ℓ₀-GIF	32.1492	0.0255	0.9751
	MM-SA-GIF	33.6431	0.00215	0.9602
[0, 150°]	TV	35.3146	0.0177	0.9830
	ℓ₀	35.5243	0.0173	0.9804
	ℓ₀-GIF	34.3459	0.0198	0.9809
	MM-SA-GIF	35.7274	0.0169	0.9617

Figure 5 shows the reconstruction results of TV minimization, LO-TV method, and MM-SA-GIF method. The second column displays the reconstruction results using the LO-TV method. Due to the lexicographic optimization method requiring optimization according to target priorities, the SART algorithm is first used for optimization. The reconstructed image obtained by SART contains numerous stripe artifacts. Subsequently, TV regularization is applied based on the SART reconstruction results, which helps to reduce some artifacts and achieve a smoother appearance. However, it is evident from the image that this method may lead to oversmooth, resulting in the loss of fine details.

Fig. 5

The reconstruction results of simulated NCAT experiment of TV minimization, LO-TV method, and MM-SA-GIF method. From the first row to the second row correspond to scanning angles of [0, 120°] and [0, 150°], respectively. From the first column to the third column, the reconstruction results of TV minimization, LO-TV method, and MM-SA-GIFmethod are in turn.The display window for the reconstruction images is [0,1]. The rectangles are the enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored.

3.2 Edam cheese experimental results

The real projection data of a sufficiently thin slice of Edam cheese [39], which has been carved with CT letters, are used to verify the superiority of the proposed MM-SA-GIF algorithm. The geometrical scanning parameters are shown in Table 4.

Table 4
The geometrical scanning parameters of CT system

Distance from source to rotation center 404.3 mm

Distance from detector to rotation center 143.5 mm

Total number of detectors 2240

The interval of projection view 1°

Pixel size of reconstruction images 0.322894 mm

Size of reconstructed image 256×256

In the experiments, the scanning ranges [0, 90°], [0, 120°] and [0, 150°] are investigated, and the reconstruction parameters are shown as follows: For TV minimization: the parameters N_TV = 20, α_TV = 0.2; For ℓ₀ minimization: (1) for the scanning range [0, 90°], the smoothing parameter η = 0.0000085, κ = 15; (2) for the scanning range [0, 120°], η = 0.0000085, κ = 15; (3) for the scanning range is [0, 150°], η = 0.0000065, κ = 15. For ℓ₀-GIF method: (1) for the scanning range [0, 90°], η = 0.0000035, κ = 15, r = 2, ɛ = 1.1 × 10^-6; (2) for the scanning range [0, 120°], η = 0.0000065, κ = 15, r = 2, ɛ = 1.1 × 10^-6; (3) for the scanning range [0, 150°], η = 0.0000055, κ = 15, r = 2, ɛ = 1.1 × 10^-5. For the above three methods, the stopping condition is reaching the number of iterations: N = 7000 for the scanning ranges [0, 90°], N = 7000 for [0, 120°], and N = 5000 for [0, 150°]. For LO-TV method:(1) for the scanning range [0, 90°], the N_SART = 500, N_TV = 350, α_TV = 0.002; (2) for the scanning range [0, 120°], the N_SART = 500, N_TV = 350, α_TV = 0.0009; (3) for the scanning range [0, 150°], the N_SART = 500, N_TV = 350, α_TV = 0.0009; For our method MM-SA-GIF, the reconstruction parameters are given in the following Table 5.

Table 5

The reconstruction parameters of MM-SA-GIF of Edam cheese Experiment

Scanning range	Parameters
[0, 90°]	J = 1400, K₁ = 5, K₂ = 100, η = 0.00000095, κ = 15, r = 2, ɛ = 1.1 × 10^-6, T₁ = 50, T₂ = 50
[0, 120°]	J = 700, K₁ = 10, K₂ = 200, η = 0.0000095, κ = 15, r = 2, ɛ = 1.1 × 10^-6, T₁ = 20, T₂ = 20
[0, 150°]	J = 500, K₁ = 10, K₂ = 100, η = 0.0000055, κ = 15, r = 2, ɛ = 1.1 × 10^-5, T₁ = 20, T₂ = 20

Figure 6 shows the reference image, which we reconstructed the image from scanning range [0, 360°] by SART method. Figures 7–9 show the results of reconstructed images by TV minimization, ℓ₀ minimization, ℓ₀-GIF method and our method for the scanning range [0, 90°], [0, 120°], [0, 150°] respectively. As can shown in Fig. 7, for the limited-angle scanning range of [0, 90°], our proposed method is clearly better at suppressing artifacts and preserving edge structures compared to TV minimization, ℓ₀ minimization, and ℓ₀-GIF method. From Fig. 8, with the increase in scanning angle, for the limited-angle scanning range is [0, 120°], the performance of all four methods in restoring image structure and suppressing artifacts are improved, but the TV minimization and ℓ₀ minimization still have obvious artifacts at the edges, and the image reconstructed by TV minimization appeared oversmoothing. The images reconstructed by ℓ₀-GIF method and our proposed method are superior to the previous two methods in artifact suppression and edge structure restoration, and from the enlarged images that our proposed method is superior to ℓ₀-GIF method in edge detail restoration. As can be seen from Fig. 9, for the scanning range is [0, 150°], due to sufficient projection data has been collected, these four image reconstruction methods can completely reconstruct images and suppress artifacts. As can be seen from the enlarged images of each picture in Fig. 9, the edge details of images reconstructed by ℓ₀ minimization are not as good as the other three methods, and compared with the other three methods, our proposed method has slight advantages in restoring image edge details.

Fig. 6

The reference image of Edam cheese.

Fig. 7

The results of Edam cheese experiment reconstructed images using different methods for the scanning range [0, 90°]. The top left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the top left and bottom right corners of each image represent the regions of interest. The red arrows indicate the restoration of edge details in the reconstructed images. The display window for the reconstruction images is [0,0.05].

Fig. 8

The results of Edam cheese experiment reconstructed images using different methods for the scanning range [0, 120°]. The top left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method.The rectangles in the top left and bottom right corners of each image represent the region of interest.The red arrows indicate whether the edge details of the reconstructed image are restored. The display window for reconstruction image is [0,0.05].

Fig. 9

The results of Edam cheese experiment reconstructed images using different methods for the scanning range [0, 150°]. The top left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the top left and bottom right corners of each image represent the region of interest.The red arrows indicate whether the edge details of the reconstructed image are restored. The display window for reconstruction image is [0,0.05].

Figure 10 shows the horizontal profiles (100th row) of the reconstructed image of Edam cheese experiment by different algorithms from different scanning ranges [0, 90°], [0, 120°], [0, 150°]. As see from Fig. 10, the image curve of our proposed method is closer to the reference image curve. The results demonstrate that our method is superior in reconstructing image details and preserving image information, indicating that the image reconstructed by our method is closer to the real image.

Fig. 10

Horizontal profiles (100th row) of the reconstructed images of Edam cheese experiment by different algorithms in different scanning ranges [0, 90°], [0, 120°], [0, 150°].

Table 6 shows the PSNR, RMSE and SSIM values of reconstruction image. As can be seen from the values, except for ℓ₀ minimization owns the highest SSIM values for scanning range [0, 150°], the proposed MM-SA-GIF method outperforms the other three methods in these indices. Additionally, the rectangles in the left top of the reconstruction image as region of interest (ROI) 1 and the rectangles in the right bottom as ROI 2. To evaluate the restoration of image edges more intuitively, we calculated the image quality index of the ROIs in the reconstructed images.

Table 6

The image quality index values of Edam cheese experiment reconstruction iamge

Scanning range	Algorithm	PSNR	RMSE	SSIM
[0, 90°]	TV	21.7292	0.0053	0.9673
	ℓ₀	23.7264	0.0042	0.9806
	ℓ₀-GIF	24.0283	0.0041	0.9801
	MM-SA-GIF	25.4876	0.0035	0.9868
[0, 120°]	TV	25.1586	0.0036	0.9829
	ℓ₀	24.5104	0.0039	0.9831
	ℓ₀-GIF	25.4682	0.0035	0.9854
	MM-SA-GIF	26.5217	0.0031	0.9887
[0, 150°]	TV	27.2508	0.0028	0.9900
	ℓ₀	27.2231	0.0028	0.9902
	ℓ₀-GIF	27.0658	0.0029	0.9892
	MM-SA-GIF	27.3401	0.0028	0.9901

Table 7 shows the image quality index values of the ROI of the reconstruction image using different methods. As seen from Table 7, the MM-SA-GIF method owns the best PSNR, RMSE and SSIM values of ROI 1 for all three scanning ranges. For ROI2, the MM-SA-GIF method owns the best SSIM value for scanning range [0, 90°] and [0, 120°]. The ℓ₀ minimization owns the best PSNR and RMSE values of ROI 2 for scanning range [0, 150°], the MM-SA-GIF method owns the best SSIM value for scanning range [0, 150°]. From the perspective of image quality indices, although the ℓ₀ minimization owns the best PSNR and RMSE values for scanning range [0, 150°], the reconstructed image of the proposed MM-SA-GIF method is more similar to the reference image in structure due to we introduce SSIM to guide the optimization.

Table 7

The image quality index values of ROI

	Scanning range		TV	ℓ₀	ℓ₀-GIF	MM-SA-GIF
ROI 1	[0, 90°]	PSNR	15.0037	15.8614	15.0135	19.7360
		RMSE	0.1778	0.1610	0.1776	0.1031
		SSIM	0.4735	0.6375	0.5818	0.7334
	[0, 120°]	PSNR	16.5481	14.8100	16.4378	19.9304
		RMSE	0.1488	0.1818	0.1507	0.1008
		SSIM	0.5284	0.5751	0.6224	0.6622
	[0, 150°]	RSNR	17.6738	19.3515	18.5944	19.5665
		RMSE	0.1307	0.1078	0.1176	0.1051
		SSIM	0.6281	0.6655	0.5998	0.6958
ROI 2	[0, 90°]	PSNR	17.0494	16.9724	20.3800	22.8729
		RMSE	0.1405	0.1417	0.0957	0.0718
		SSIM	0.4763	0.5652	0.6156	0.6614
	[0, 120°]	PSNR	16.7780	18.0912	19.2743	21.7916
		RMSE	0.1449	0.1246	0.1087	0.0814
		SSIM	0.5175	0.5784	0.5876	0.6056
	[0, 150°]	PSNR	17.2231	25.0766	19.7024	20.9076
		RMSE	0.1377	0.0557	0.1035	0.0901
		SSIM	0.6169	0.6698	0.6433	0.7112

Figure 11 presents the reconstruction results of the TV minimization, LO-TV method, and MM-SA-GIF method. The second column showcases the reconstruction results using the LO-TV method. It is evident from the image that the reconstructed image obtained through the LO-TV algorithm exhibits fold artifacts. This occurs because the SART, when used with insufficient projection data, often produces significant artifacts. Consequently, applying TV regularization on top of the SART reconstruction results can still result in fold artifacts in the reconstructed image. In addition, this method often leads to oversmoothing, which causes the loss of details and reduces the overall clarity of the image.

Fig. 11

The reconstruction results of Edam cheese experiment of TV minimization, LO-TV method, and MM-SA-GIF method. From the first row to the third row correspond to scanning angles of [0, 90°], [0, 120°] and [0, 150°], respectively. From the first column to the third column, the reconstruction results of TV minimization, LO-TV method, and MM-SA-GIFmethod are in turn.The display window for the reconstruction images is [0,0.05]. The rectangles are the enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored.

3.3 Drumstick experimental results

Another real projection data of a drumstick is utilized to evaluate the perfomance of the proposed method. The geometrical scanning parameters of CT systems are shown in Table 8.

Table 8
The geometrical scanning parameters of CT system

Distance from source to rotation center 315.7 mm

Distance from detector to rotation center 879.3 mm

Total number of detectors 2048

The interval of projection view 0 . 3519°

Pixel size of reconstruction images 0.2112 mm

Size of reconstructed image 512×512

In the experiments, we investigated scanning ranges of [0, 90°], [0, 120°], and [0, 150°], with the following reconstruction parameters: For TV minimization: the parametes N_TV = 20, α_TV = 0.2. For ℓ₀ minimization: (1) for the scanning range [0, 90°], η = 0.00000085, κ = 15; (2) for the scanning range [0, 120°], η = 0.00000085, κ = 15; (3) for the scanning range is [0, 150°], η = 0.00000075, κ = 15. For ℓ₀-GIF method: (1) for the scanning range [0, 90°], η = 0.000005, κ = 15, r = 2, ɛ = 1.1 × 10^-6; (2) for the scanning range [0, 120°], η = 0.00000085, κ = 15, r = 2, ɛ = 1.1 × 10^-6; (3) for the scanning range [0, 150°], η = 0.00000075, κ = 15, r = 2, ɛ = 1.1 × 10^-6. The stopping condition for the above three methods is reaching the number of iterations N = 3000 for each respective scanning range. For LO-TV method:(1) for the scanning range [0, 90°], the N_SART = 500, N_TV = 350, α_TV = 0.0005; (2) for the scanning range [0, 120°], the N_SART = 500, N_TV = 350, α_TV = 0.0005; (3) for the scanning range [0, 150°], the N_SART = 500, N_TV = 350, α_TV = 0.0004; For our method MM-SA-GIF, the reconstruction parameters are provided in Table 9.

Table 9

The reconstruction parameters of Drumstick Experiment

Scanning range	Parameters
[0, 90°]	J = 150, K₁ = 10, K₂ = 50, η = 0.0000002, κ = 15, r = 2, ɛ = 1.1 × 10^-6, T₁ = 20, T₂ = 20
[0, 120°]	J = 100, K₁ = 10, K₂ = 100, η = 0.00000045, κ = 15, r = 2, ɛ = 1.1 × 10^-6, T₁ = 20, T₂ = 20
[0, 150°]	J = 100, K₁ = 10, K₂ = 50, η = 0.00000025, κ = 15, r = 2, ɛ = 1.1 × 10^-6, T₁ = 20, T₂ = 20

Figure 12 shows the reference image, which we reconstructed the image by TV minimization from the scanning range [0, 360°]. Figures 14–16 display the reconstructed drumstick images obtained using TV minimization, ℓ₀ minimization, ℓ₀-GIF method and our method for scanning range [0, 90°], [0, 120°] and [0, 150°]. As can be seen from Fig. 14, for the scanning range [0, 90°], our proposed method outperforms TV minimization, ℓ₀ minimization and ℓ₀-GIF method in suppressing artifacts and preserving edges, but the image reconstructed by our proposed method is oversmooth. From Fig. 15, for the scanning range is [0, 120°], the reconstructed image of TV minimization and ℓ₀ minimization have streak artifacts and the edges are not clear enough, with ℓ₀-GIF method showing better artifact suppression compared to ℓ₀ minimization.In the enlarged images at the bottom right, our proposed method effectively suppresses artifacts while preserving image edges, particularly evident in soft tissue edges. Moreover, in the enlarged images at the top left, our proposed method demonstrates superior artifact suppression in detail compared to other methods. As can be seen from Fig. 16, for the scanning range is [0, 150°], the reconstructed image of ℓ₀ minimization is still exists a few streak artifacts, and the edge details are not clear enough. Compared with ℓ₀ minimization, the other three method can effectively suppress artifacts. From the enlarged detail images in the right bottom, our proposed method is better in restoring image edge details.

Fig. 12

The reference image of drumstick.

Fig. 13

Horizontal profiles (256th row) of the reconstructed images of drumstick experiment by different algorithms in different scanning ranges [0, 90°], [0, 120°], [0, 150°]. The rectangles in the left up are the details of the horizontal profiles by red arrows.

Fig. 14

The results of drumstick experiment reconstructed images using different methods for the scanning range [0, 90°]. The top Left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the bottom right corners are the enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored.

Fig. 15

The results of drumstick experiment reconstructed images using different methods for the scanning range [0, 120°]. The top Left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the top left and bottom right corners are the enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored.

Fig. 16

The results of drumstick experiment reconstructed images using different methods for the scanning range [0, 150°]. The top Left is the result reconstructed by TV minimization, the top right image using ℓ₀ minimization, the bottom left image using ℓ₀-GIF method, and the bottom right image using the MM-SA-GIF method. The rectangles in the bottom right corners are the enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored.

Figure 13 displays the horizontal profiles (256th row) of the reconstructed image of drumstick experiment by different algorithms from different scanning ranges [0, 90°], [0, 120°], [0, 150°]. The rectangles in the left up are the enlarged image of horizontal profiles. As can be seen from Fig. 13, for the scanning ranges are [0, 90°] and [0, 120°], the horizontal profiles of the image of our proposed method is closer to the horizontal profiles of the reference image. For the scanning range is [0, 150°], the horizontal profiles of these four methods are all close to the reference image, but as can be seen from the enlarged image, our proposed method are more close to that of the reference image. It means that the image reconstructed by our proposed method is closer to the reference image in structure.

Table 10 shows the PSNR, RMSE and SSIM values of reconstructed images by different methods. In Table 10, for the scanning ranges of [0, 90°] and [0, 120°], the PSNR, RMSE and SSIM values of the proposed MM-SA-GIF method are better than those of the other three methods. For the scanning range of [0, 150°], the PSNR, RMSE and SSIM values of the TV minimization are better than other methods, however, as can be seen from the Fig. 16, the proposed method MM-SA-GIF is superior to TV minimization in suppressing artifacts, preserving reconstruction image edged and restoring image structural details. We should note that the values of the image quality index will be affected by the reference image, and because we introduce the SA algorithm to accept the degradation solution in the form of probability, so the SSIM values of TV minimization may be better than our proposed method. However, the SSIM values will not decline too much, because the probability of acceptance will be very low if the SSIM values decline too much.

Table 10

The image quality index values of drumstick experiment reconstruction iamge

Scanning range	Algorithm	PSNR	RMSE	SSIM
[0, 90°]	TV	25.5208	0.0035	0.9758
	ℓ₀	25.3099	0.0036	0.9780
	ℓ₀-GIF	25.2980	0.0036	0.9746
	MM-SA-GIF	25.7064	0.0034	0.9793
[0, 120°]	TV	26.9861	0.0030	0.9828
	ℓ₀	26.7298	0.0031	0.9834
	ℓ₀-GIF	26.8199	0.0030	0.9825
	MM-SA-GIF	27.0858	0.0029	0.9843
[0, 150°]	TV	27.8255	0.0027	0.9860
	ℓ₀	27.3350	0.0028	0.9856
	ℓ₀-GIF	27.4549	0.0028	0.9853
	MM-SA-GIF	27.5180	0.0028	0.9858

Figure 17 shows the reconstruction results of the TV minimization, LO-TV method, and MM-SA-GIF method. The second column displays the reconstruction results using the LO-TV method. Similar to the results from the simulated NCAT experiment and the Edam-cheese experiment, the images obtained using the LO-TV method exhibit significant stripe artifacts and continue to suffer from oversmoothing.

Fig. 17

The reconstruction results of drumstick experiment of TV minimization, LO-TV method, and MM-SA-GIF method. From the first row to the third row correspond to scanning angles of [0, 90°], [0, 120°] and [0, 150°], respectively. From the first column to the third column, the reconstruction results of TV minimization, LO-TV method, and MM-SA-GIFmethod are in turn. The rectangles are the enlarged detail areas. The red arrows indicate whether the edge details of the reconstructed image are restored.

4 Conclusions and perspectives

In this paper, we proposed a multiobjective optimization model for limited-angle CT reconstruction, aiming to enhance image quality. The multiobjective optimization model incorporates two objective functions: data fidelity term and ℓ₀-norm of the image gradient. Unlike traditional scalarzation medthods minimizing the objective functions, we use the image quality index SSIM to guide the optimization process and solve this model by an iterative method. Firstly, we use SART to optimize the data fidelity term, and introduce SSIM to guided the optimization process. We cancluate the SSIM values between the iterative points and reference image. When the image quality degrades, the SA algorithm is introduced to accept the degradation solution in the form of probability. If the degraded solution is accepted, we repeat the use of SART to optimize the data fidelity term. Otherwise, we introduce the GIF method to further preserve the image edges and generate alternative solutions. Secondly, we take the result of the first step as the initial point and add it as a constraint to the second objective function. We use the ℓ₀ minimization to optimize ℓ₀-norm of the image gradient, and introduce SSIM, the SA algorithm, and GIF to guide the optimization process, similar to the first step. In the experiment part, we validate the superiority of our proposed method using both simulated data and real projection data experiments. The reconstructed results demonstrate that our method outperforms other methods in preserving image edges and suppressing artifacts. Horizontal profiles of the reconstructed images by different methods in various scanning ranges are plotted to visually observe the differences between the reconstructed image and the reference image. Compared with other methods, the profiles of our method are closer to those of the reference image than those of the other methods. Additionally, it is important to note that the choice of reference image can significantly impact the experimental data, and thus, careful selection of the reference image is necessary.

While our proposed method has advantages in suppressing artifacts and preserving edge details, the parameters selected in the reconstruction process rely on experience for manual debugging, which is laborious and time-consuming. If the scanning range is too small, the reconstructed image by our method maybe oversmooth. Furthermore, we need to acknowledge that our algorithm may have a lower SSIM value than other algorithms. This is primarily due to our utilization of the SA algorithm, which probabilistically accepts degraded solutions, leading to potentially lower SSIM values. However, despite this, the reconstructed images produced by our algorithm significantly outperform other algorithms in terms of edge preservation and artifact suppression.

In future research, since the SSIM value is influenced by the reference image and thereby impacts our reconstruction process. In future work, we plan to conduct a more comprehensive analysis to further understand the impact of different reference images on the SSIM values. Additionally, we will explore whether there is an image quality index that does not rely on a reference image.

Footnotes

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NSFC) (No.61801086, No.62371184, No.11991024, No.12171063 and No.11671062); Natural Science Foundation of Chongqing (No.cstc2019jcyj-msxmX0345, No.CSTB2023NSCQ-MSX0148, and No.cstc2021jcyj-msxmX0280); The Project for Creative Research Groups at Institutions of Higher Education in Chongqing (No.CXQT20014). Scientific and technological research program of Chongqing Municipal Education Commission (KJQN202300535, KJQN202000534, KJZD-K202100503, KJZD-K202300506 and KJQN202100521).

References

Tingberg

, X-ray tomosynthesis: a review of its use for breast andchest imaging, Radiation Protection Dosimetry 139(1–3) (2010), 100–107. doi:https://doi.org/10.1093/rpd/ncq099

Bachar

, Siewerdsen

, Daly

, Jaffray

, Irish

, Image quality and localization accuracy in c-arm tomosynthesis-guided head and neck surgery, Medical Physics 34(12) (2007), 4664–4677. doi:https://doi.org/10.1118/1.2799492

Cho

, Kim

, Youn

, Kim

, A feasibility study of digitaltomosynthesis for volumetric dental imaging, Journal of Instrumentation 7(03) (2012), P03007. doi:https://doi.org/10.1088/1748-0221/7/03/P03007

Natterer

, Natterer

, Computerized tomography, The Mathematics of Computerized Tomography (1986), 1–8.

Willemink

M.J.

, Noël

P.B.

, The evolution of image reconstruction for ct—from filtered back projection to artificial intelligence, European Radiology 29 (2019), 2185–2195. doi:https://doi.org/10.1007/s00330-018-5810-7

Beister

, Kolditz

, Kalender

W.A.

, Iterative reconstruction methods in x-ray ct, Physica Medica 28(2) (2012), 94–108. doi:https://doi.org/10.1016/j.ejmp.2012.01.003

Andersen

A.H.

, Kak

A.C.

, Simultaneous algebraic reconstruction technique (sart): a superior implementation of the art algorithm, Ultrasonic Imaging 6(1) (1984), 81–94. doi:https://doi.org/10.1177/016173468400600107

Dempster

A.P.

, Laird

N.M.

, Rubin

D.B.

, Maximum likelihood from incomplete data via the em algorithm, Journal of the Royal Statistical Society: Series B (methodological) 39(1) (1977), 1–22. doi:https://doi.org/10.1111/j.2517-6161.1977.tb01600.x

Sidky

E.Y.

, Kao

C.-M.

, Pan

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam ct, Journal of X-ray Science and Technology 14(2) (2006), 119–139. doi:https://doi.org/10.48550/arXiv.0904.4495

10.

Sidky

E.Y.

, Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total variation minimization, Physics in Medicine & Biology 53(17) (2008), 4777. doi:https://doi.org/10.1088/0031-9155/53/17/021

11.

Chen

, Jin

, Li

, Wang

, A limited-angle ctreconstruction method based on anisotropic tv minimization, Physics in Medicine & Biology 58(7) (2013), 2119. doi:https://doi.org/10.1088/0031-9155/58/7/2119

12.

Gong

, Zeng

, Wang

, Image reconstruction model forl imited-angle ct based on prior image induced relative total variation, Applied Mathematical Modelling 74 (2019), 586–605. doi:https://doi.org/10.1016/j.apm.2019.05.020

13.

Luo

, Li

, Tu

, Wu

, Adaptive weighted total variation minimization based alternating direction method of multipliers for limited angle ct reconstruction, IEEE Access 6 (2018), 64225–64236. doi:https://doi.org/10.1109/ACCESS.2018.2873713

14.

Liu

, Ma

, Fan

, Liang

, Adaptive-weighted total variation minimization for sparse data toward lowdose x-ray computed tomography image reconstruction, Physics in Medicine & Biology 57(23) (2012), 7923. doi:https://doi.org/10.1088/0031-9155/57/23/7923

15.

Wang

, Nakamoto

, Zhang

, Liu

, Reweighted anisotropic total variation minimization for limited angle ct reconstruction, IEEE Transactions on Nuclear Science 64(10) (2017), 2742–2760. doi:https://doi.org/10.1109/TNS.2017.2750199

16.

Xiao

, Zhong

, Li

, Hu

, Xu

, A new limited-angle ct reconstruction algorithm based on the local anisotropic totalvariation restoration of continuity, Journal of Instrumentation 17(12) (2022), P12018. doi:https://doi.org/10.1088/1748-0221/17/12/P12018

17.

, Zeng

, ℓ₀ gradient minimization based imagereconstruction for limited-angle computed tomography, PloS One 10(7) (2015), e0130793. doi:https://doi.org/10.1371/journal.pone.0130793

18.

Dong

, Zhang

, An efficient algorithm for .0 minimization inwavelet frame based image restoration, Journal of Scientific Computing 54 (2013), 350–368. doi:https://doi.org/10.1007/s10915-012-9597-4

19.

Zhang

, Dong

, Lu

, ℓ₀ minimization for wavelet framebased image restoration, Mathematics of Computation 82(282) (2013), 995–1015. doi:https://doi.org/10.1090/s0025-5718-2012-02631-7

20.

Zhou

, Cai

J.-F.

, Gao

, Adaptive tight frame based medicalimage reconstruction: a proof-of-concept study for computed tomography, Inverse Problems 29(12) (2013), 125006. doi:https://doi.org/10.1088/0266-5611/29/12/125006

21.

Wang

, Zeng

, Guo

, Zhang

, Wavelet tight frame and prior image-based image reconstruction from limited-angle projection data, Inverse Problems & Imaging 11(6) (2017), 917–948. doi:https://doi.org/10.3934/ipi.2017043

22.

Zhao

, Gao

, Ding

, Molloi

, Tight-frame based iterative image reconstruction for spectral breast ct, Medical Physics 40(3) (2013), 031905. doi:https://doi.org/10.1118/1.4790468

23.

, Sun

, Tang

, Guided image filtering, IEEE Transactions on Pattern Analysis and Machine Intelligence 35(6) (2013), 1397–1409. doi:https://doi.org/10.1109/TPAMI.2012.213

24.

, Qu

, Liu

, Simultaneous algebraic reconstruction technique based on guided image filtering, Optics Express 24(14) (2016), 15897–15911. doi:https://doi.org/10.1364/OE.24.015897

25.

Wang

, Wang

, Zhao

, Yu

, Huang

, Guided image filtering based ℓ₀ gradient minimization for limited-angle ct image reconstruction, Journal of Inverse and Ill-posed Problems 29(4) (2021), 587–598. doi:https://doi.org/10.1515/jiip-2020-0096

26.

Wang

, Yue

, Wang

, Yu

, Image reconstruction method for limited-angle ct based on total variation minimization using guided image filtering, Medical & Biological Engineering & Computing 60(7) (2022), 2109–2118. doi:https://doi.org/10.1007/s11517-022-02579-z

27.

Han

, Bian

, Ritman

E.L.

, Sidky

E.Y.

, Pan

, Optimization-based reconstruction of sparse images from few-view projections, Physics in Medicine & Biology 57(16) (2012), 5245. doi:https://doi.org/10.1088/0031-9155/57/16/5245

28.

Han

, Bian

, Eaker

D.R.

, Kline

T.L.

, Sidky

E.Y.

, Ritman

E.L.

, Pan

, Algorithm-enabled low-dose microct imaging, IEEE Transactions on Medical Imaging 30(3) (2010), 606–620. doi:https://doi.org/10.1109/TMI.2010.2089695

29.

Zhang

, Chen

, Xia

, Sidky

E.Y.

, Pan

, Directional-tv algorithm for image reconstruction from limited angular-range data, Medical Image Analysis 70 (2021), 102030. doi:https://doi.org/10.1016/j.media.2021.102030

30.

Sawaragi

, Nakayama

, Tanino

, Theory of multiobjective optimization, Elsevier, 1985.

31.

Ehrgott

, Multicriteria optimization, Vol. 491, Springer Science& Business Media, 2005.

32.

Miettinen

, Some methods for nonlinear multi-objective optimization, in: International conference on evolutionary multi-criterion optimization, Springer, 2001, pp. 1.20. doi:https://doi.org/10.1007/3-540-44719-9_1

33.

Wang

, A perspective on deep imaging, IEEE Access 4 (2016), 8914–8924. doi:https://doi.org/10.1109/ACCESS.2016.2624938

34.

Fan

, Shan

, Kalra

M.K.

, Singh

, Qian

, Getzin

, Teng

, Hahn

and Wang

, Quadratic autoencoder (q-ae) for low-dose ct denoising, IEEE Transactions on Medical Imaging 39(6) (2019), 2035–2050. doi:https://doi.org/10.1109/TMI.2019.2963248

35.

Chen

, Zhang

, Kalra

M.K.

, Lin

, Chen

, Liao

, Zhou

, Wang

, Low-dose ct with a residual encoder-decoder convolutional neural network, IEEE Transactions on Medical Imaging 36(12) (2017), 2524–2535. doi:https://doi.org/10.1109/TMI.2017.2715284

36.

Yang

, Yan

, Zhang

, Yu

, Shi

, Mou

, Kalra

M.K.

, Zhang

, Sun

and Wang

, Low-dose ct image denoising using agenerative adversarial network with wasserstein distance andperceptual loss, IEEE Transactions on Medical Imaging 37(6) (2018), 1348–1357. doi:https://doi.org/10.1109/TMI.2018.2827462

37.

Sidky

E.Y.

, Lorente

, Brankov

J.G.

, Pan

, Do cnns solve thect inverse problem? IEEE Transactions on Biomedical Engineering 68(6) (2021), 1799–1810. doi:https://doi.org/10.1109/TBME.2020.3020741

38.

Kirkpatrick

, Gelatt

C.D.

Jr , Vecchi

M.P.

, Optimization by simulated annealing, Science 220(4598) (1983), 671–680. doi:https://doi.org/10.1126/science.220.4598.671

39.

Bubba

T.A.

, Juvonen

, Lehtonen

, März

, Meaney

, Purisha

and Siltanen

, Tomographic x-ray data of carved cheese, arXiv preprint arXiv:1705.05732 (2017). doi:https://doi.org/10.48550/arXiv.1705.05732

Distance from source to rotation center	404.3 mm
Distance from detector to rotation center	143.5 mm
Total number of detectors	2240
The interval of projection view	1°
Pixel size of reconstruction images	0.322894 mm
Size of reconstructed image	256×256

Distance from source to rotation center	315.7 mm
Distance from detector to rotation center	879.3 mm
Total number of detectors	2048
The interval of projection view	0 . 3519°
Pixel size of reconstruction images	0.2112 mm
Size of reconstructed image	512×512