A directional relative TV algorithm for sparse-view CT reconstruction

Abstract

Objective:

Computed tomography (CT) is a widely used medical imaging modality, but its radiation exposure poses potential risks to human health. Sparse-view scanning has emerged as an effective approach to reduce radiation dose; however, images reconstructed using the filtered back-projection (FBP) algorithm from sparse-view projections often suffer from severe streak artifacts. Achieving high-quality CT image reconstructed from sparse-view projections remains a challenging task.

Methods:

Building on compressed sensing (CS), the total variation (TV) algorithm is applied for high-quality sparse-view reconstruction. We further propose a relative total variation (RTV) algorithm to enhance the accuracy of sparse-view reconstruction. Experimental results indicate that while the RTV algorithm improves accuracy, it has limitations in edge preservation. To address this, inspired by the success of directional TV (DTV) in limited-angle reconstruction, we develop a directional relative TV (DRTV) model. This model applies the RTV technique in both x and y directions independently, and we derive its adaptive steepest descent projection onto convex set (ASD-POCS) solution algorithm.

Results:

Experiments conducted on simulated phantoms and real CT images demonstrate the correctness, convergence, and superior performance of the DRTV algorithm in sparse-view reconstruction. Compared with the TV, DTV, and RTV algorithm, the DRTV algorithm exhibits superior preservation of structural features and texture details.

Significance:

The DRTV algorithm represents an advanced method for high-precision sparse-view CT reconstruction, providing stable and accurate results. Moreover, the approach is applicable to other medical imaging modalities.

Keywords

computed tomography sparse-view reconstruction relative TV directional TV adaptive steepest descent projection onto convex set algorithm

Introduction

With the advantages of high resolution, fast imaging speed and versatility, computed tomography (CT) has become a widely used imaging modality in clinical diagnosis and treatment nowadays. However, the health risks of radiation exposure from CT should not be underestimated. Sparse-view scanning is a configuration that effectively reduces CT radiation dose. But images reconstructed using the filtered back-projection (FBP) algorithm,¹ a traditional and commercially available reconstruction method, from sparse-view projections often suffer from severe streak artifacts.² How to achieve high-precision sparse-view CT reconstruction is a challenging problem.

With the aid of compressive sensing (CS) theory,³ many ill-posed inverse problems, including sparse-view reconstruction, can now be effectively solved. Among the most representative solutions is the development of total variation (TV) type algorithms. In 2006, the model problem of reconstructing an object from incomplete frequency samples was studied by Candes et al., where it was shown that function can be reconstructed exactly as the solution to the minimization problem.⁴ That is to say that exact recovery may be obtained by solving a convex optimization problem. On this basis, Sidky et al. develop and investigate an iterative image reconstruction algorithm based on the minimization of the image TV that applies to divergent-beam CT, which can realize sparse reconstruction with high accuracy.⁵ Subsequently, in 2008, Sidky et al. further optimized the algorithm and proposed the adaptive steepest descent projection onto convex set (ASD-POCS) algorithm.⁶

Several variants of TV have since emerged to address the shortcomings of TV algorithms for specific problems. In 2010, Synho Do et al. present an improved regularization technique by incorporating higher-order (HOTV) derivatives to reduce staircase artifacts without sacrificing edge sharpness.⁷ Later people also applied HOTV in other scenarios.^8,9 In 2011, since the low-contrast structures tend to be smoothed out by the TV regularization, therefore Tian et al. developed an iterative CT reconstruction algorithm with edge-preserving TV (EPTV) regularization to reconstruct CT images from highly under-sampled data obtained at low mAs levels.¹⁰ In 2012, since a conventional TV minimization algorithm often suffers from over-smoothness on the edges of the resulting image, Liu et al. proposed an adaptive-weighted TV (AwTV) minimization algorithm derived by considering the anisotropic edge property among neighboring image voxels.¹¹ In 2012, Xu et al. proposed the relative total variation (RTV) model, which can extract meaningful structures under the complication of texture patterns.¹² In 2013, since the TV minimization process is isotropic, suggesting that it is unfit for limited-angle CT. Therefore Chen et al. proposed an anisotropic TV (ATV) minimization method.¹³ In the same year, Ning et al. proposed a new image reconstruction algorithm based on lp norm compressive sensing by combining the penalty function and revised Hesse sequence quadratic programming, and using block compressive sensing.¹⁴ In 2014, Sidky et al. proposed the constrained total p variation (TpV) model and solved the model using the Chambolle-Pock (CP)^15–17 optimization algorithm.¹⁸ In 2015, Rigie et al. proposed the total nuclear variation (TnV) model, which better preserves the boundary information by encouraging the different channels to have the same edge structure and the gradient vectors to point to the same direction.¹⁹ In 2017, Wang et al. proposed a new iteratively reweighted ATV method, in which a reweighted technique is incorporated into the idea of ATV.²⁰ In 2021, Zhang et al. proposed the directional total variation (DTV) model, which introduces the direction of the gradient to constrain the variation of the image, thus preserving the structural features and texture details of the image.²¹ In the same year, in ASD-POCS algorithm, the existing gradient expression of the TV-type norm appears too complicated in the implementation code and reduces image reconstruction speed. To address this issue, Qiao developed a simple and fast ASD-POCS algorithm.²² In the same year, Qiao et al. proposed the balanced total variation (bTV) model, which can guarantee convergence and achieve fast convergence in 3D electron paramagnetic resonance imaging (EPRI)^23–25 applications.²⁶ In 2023, Qiao et al. proposed a data divergence constrained, total nuclear variation minimization model and its CP solving algorithm.²⁷ In 2024, Liu et al. proposed an edge-preserving total nuclear variation (EPTV_N) minimization algorithm for sparse reconstruction in EPRI.²⁸ In 2024, Qiao et al. proposed to use the DTV algorithm in sparse reconstruction for EPR imaging, which outperforms existing FBP and TV-like algorithms as well as deep learning-based methods.²⁹ In the same year, Fang et al. used the DTV algorithm for fast EPR imaging, which achieves a 10-fold speedup compared to the standard FBP algorithm.³⁰

The TV algorithm is not conducive to preserving the image structure during image reconstruction as it penalizes all the gradients in the reconstructed image and may also produce streak artifacts. The RTV model can be used to penalize different sizes of image gradients, which can adaptively protect the image structure in sparse projection CT reconstruction and help to suppress the blocky artifacts to a certain extent.^31–33 The DTV model imposes directional TV constraints on the image and accurately recovers the phantoms from data generated over a significantly reduced angular range, and that it considerably diminishes artifacts observed otherwise in reconstructions of existing algorithms.^34,35 Since RTV is weak in edge-preserving ability and DTV has edge-protecting advantages, in this paper, we propose a directional relative total variation (DRTV) reconstruction model, which is a fusion of RTV and DTV rather than a simple combination of the two. It can be understood as a directional treatment of relative TV, or it can be seen as a relativized improvement of directional TV. From the point of view of sparse optimization, DRTV is able to utilize more accurate a priori information corresponding to DRTV regular term compared to RTV, thus exerting stronger constraint ability on the solution. Through this fusion, the advantages of both are fully combined, and it is expected to realize high-precision sparse reconstruction.

In summary, the main contributions of this work are listed as follows:

We combine the artifact reduction advantage of RTV and the edge protection advantage of DTV to propose a DRTV model for sparse CT reconstruction.

We derived the DRTV-ASD-POCS solution algorithm in the ASD-POCS framework.

We design simulation experiments and real CT experiments to verify the correctness, stability and superior performance of the proposed algorithm in sparse reconstruction.

The rest of the paper is organized as follows. In Section II, we design the DRTV model and derive the DRTV-ASD-POCS solving algorithm in the framework of ASD-POCS algorithm. In Section III, we organize experiments to verify the correctness and convergence of the algorithm as well as the superiority it shows on sparse reconstruction. The discussion and conclusions are given in Section IV.

Methods

Imaging system model

In this paper, two-dimensional (2D) parallel-beam CT reconstruction is used as a research object. The discrete to discrete (D2D) imaging model for 2D CT, can be expressed as:

g = A u,

(1)

where

g

is the discrete projection data, a column vector of size M, and each measurement is

g_{i}

, i=1,2,…,M; the image to be solved is represented in discrete form as

u

, a column vector of length N, and the elements in the image are

u_{j}

, j=1,2,…,N;

A

is the system matrix of size

M \times N

and the elements are denoted as

A_{i, j}

. In this paper, the system matrix is solved by Siddon's ray-driven method,³⁶ thus

A_{i, j}

denotes the length of the i th ray through the j th pixel.

Typically, the system matrix is large-scale, ill-posed and under-conditioned, so direct matrix inversion is computationally intractable. In order to improve the reconstruction accuracy, the model needs to be transformed into an optimization problem and solved using the appropriate method.

RTV model

RTV model shown as

u * = \underset{u}{argmin} ‖ u ‖_{RTV} s . t . ‖ g - A u ‖_{2} \leq ε

(2)

where

u *

denotes the optimal solution that satisfies the constraints, i.e., the reconstructed image;

‖ g - A u ‖_{2} \leq ε

is the data fidelity term, which aims to reduce the gap between the real projection

g

and the simulated projection

A u

;

ε

is the data tolerance, the size of which is related to the noise;

‖ u ‖_{RTV}

defined as follows:

‖ u ‖_{RTV} = \sum_{p} \frac{D_{x} (p)}{L_{x} (p) + ϵ} + \frac{D_{y} (p)}{L_{y} (p) + ϵ}

(3)

where the parameter

ϵ

avoids a denominator of 0;

D_{x} (p)

and

D_{y} (p)

are windowed total variation (WTV), which are defined as follows:

D_{x} (p) = \sum_{q \in R (p)} g_{p, q} \cdot | {(\partial_{x} u)}_{q} |

(4)

D_{y} (p) = \sum_{q \in R (p)} g_{p, q} \cdot | {(\partial_{y} u)}_{q} |

(5)

where

R (p)

denotes a localized rectangular region centered on pixel p, q is the pixel index within this rectangular region, and

g_{p, q}

denotes a Gaussian function with

σ

as the standard deviation, i.e.,

g_{p, q} \propto \exp (- \frac{{(x_{p} - x_{q})}^{2} + {(y_{p} - y_{q})}^{2}}{2 σ^{2}}),

(6)

where

L_{x} (p)

and

L_{y} (p)

denote windowed inherent variation (WIV), which are defined as follows:

L_{x} (p) = | \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{x} u)}_{q} |,

(7)

L_{y} (p) = | \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{y} u)}_{q} | .

(8)

DTV model

DTV model shown as

u * = \underset{u}{argmin} \frac{1}{2} ‖ g - A u ‖_{2}^{2} s . t . ‖ u ‖_{T V_{x}} \leq t_{x}, ‖ u ‖_{T V_{y}} \leq t_{y}, and u_{i} \geq 0,

(9)

where

‖ u ‖_{T V_{x}}

and

‖ u ‖_{T V_{y}}

are the DTV norms along x and y directions, respectively;

t_{x}

and

t_{y}

denote the upper bounds on them. The definitions of

‖ u ‖_{T V_{x}}

and

‖ u ‖_{T V_{y}}

are

‖ u ‖_{T V_{x}} = ‖ (| D_{x} u |) ‖_{1}, ‖ u ‖_{T V_{y}} = ‖ (| D_{y} u |) ‖_{1} .

(10)

They are the $ℓ_{1}$ norm of the directional gradient magnitude transforms along x and y directions. Matrices $D_{x}$ and $D_{y}$ of size N × N may be denoted as

(D_{x} u)_{x, y} = {\begin{matrix} u_{x, y} - u_{x - 1, y} & x \in [2, N_{x}] \\ 0 & x = 1 \end{matrix},

(11)

(D_{y} u)_{x, y} = {\begin{matrix} u_{x, y} - u_{x, y - 1} & y \in [2, N_{y}] \\ 0 & y = 1 \end{matrix} .

(12)

DRTV model

In this paper, we incorporate the idea of DTV into the RTV, which constrains the x and y directions, respectively. The proposed DRTV model defined as

u * = \underset{u}{argmin} ‖ u ‖_{DRTV} s . t . ‖ g - A u ‖_{2} \leq ε,

(13)

where the norm term

‖ u ‖_{DRTV} = ‖ u ‖_{{DRTV}_{x}} + b \cdot ‖ u ‖_{{DRTV}_{y}}

, and the parameter b is the weight to balance the intensities of the x and y directions in the reconstruction.

‖ u ‖_{{DRTV}_{x}}

and

‖ u ‖_{{DRTV}_{y}}

are defined as follows, respectively,

‖ u ‖_{{DRTV}_{x}} = \sum_{p} \frac{D_{x} (p)}{L_{x} (p) + ϵ},

(14)

‖ u ‖_{{DRTV}_{y}} = \sum_{p} \frac{D_{y} (p)}{L_{y} (p) + ϵ} .

(15)

Equations (14) and (15) can be written from equations (4), (5), (7) and (8) as

‖ u ‖_{{DRTV}_{x}} = \sum_{p} \frac{\sum_{q \in R (p)} g_{p, q} \cdot | {(\partial_{x} u)}_{q} |}{| \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{x} u)}_{q} | + ϵ},

(16)

‖ u ‖_{{DRTV}_{y}} = \sum_{p} \frac{\sum_{q \in R (p)} g_{p, q} \cdot | {(\partial_{y} u)}_{q} |}{| \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{y} u)}_{q} | + ϵ} .

(17)

By reorganizing the terms and grouping elements containing $| {(\partial_{x} u)}_{q} |$ in equation (16) shown as

\begin{aligned} ‖ u ‖_{{DRTV}_{x}} & = \sum_{q} \sum_{p \in R (q)} \frac{g_{p, q}}{| \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{x} u)}_{q} | + ϵ} | {(\partial_{x} u)}_{q} | \\ \approx \sum_{q} \sum_{p \in R (q)} \frac{g_{p, q}}{| \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{x} u)}_{q} | + ϵ} \frac{1}{| {(\partial_{x} u)}_{q} | + ϵ_{u}} (\partial_{x} u)_{q}^{2} \\ = \sum_{q} z_{x q} w_{x q} (\partial_{x} u)_{q}^{2} . \end{aligned}

(18)

The second line in equation (18) is an approximation obtained due to the numerical stability of the introduced $ϵ_{u}$ . By rearranging these terms, it is decomposed into a quadratic term $(\partial_{x} u)_{q}^{2}$ and a nonlinear component $z_{x q} w_{x q}$ which are, respectively,

z_{x q} = \sum_{p \in R (q)} \frac{g_{p, q}}{| \sum_{q \in R (p)} g_{p, q} \cdot {(\partial_{x} u)}_{q} | + ϵ} = {(G_{σ} * \frac{1}{| G_{σ} * \partial_{x} u | + ϵ})}_{q},

(19)

w_{x q} = \frac{1}{| {(\partial_{x} u)}_{q} | + ϵ_{u}} .

(20)

Equation (19) indicates that $z_{x}$ for each pixel actually incorporates neighboring gradient information in an isotropic spatial filter manner. $G_{σ}$ is a Gaussian filter with standard deviation $σ$ . $w_{x}$ is only related to the pixel-wise gradient.

Similarly, we can express equation (17) as

‖ u ‖_{{DRTV}_{y}} = \sum_{q} z_{y q} w_{y q} (\partial_{y} u)_{q}^{2},

(21)

where

(\partial_{y} u)_{q}^{2}

is the quadratic y-component partial derivative and

z_{y q} w_{y q}

similarly the non-linear part. They are respectively

z_{y q} = {(G_{σ} * \frac{1}{| G_{σ} * \partial_{y} u | + ϵ})}_{q},

(22)

w_{y q} = \frac{1}{| {(\partial_{y} u)}_{q} | + ϵ_{u}} .

(23)

Equations (14) and (15) written in matrix form can be obtained

‖ u ‖_{{DRTV}_{x}} = v_{u}^{T} C_{x}^{T} Z_{x} W_{x} C_{x} v_{u},

(24)

‖ u ‖_{DRTV y} = v_{u}^{T} C_{y}^{T} Z_{y} W_{y} C_{y} v_{u},

(25)

where

v_{u}

is the vector representation of

u

C_{x}

and

C_{y}

are the Toeplitz matrices obtained from the discrete gradient operator of the forward approximation, and

Z_{x}

Z_{y}

W_{x}

and

W_{y}

are all diagonal arrays. Their diagonal values are respectively

Z_{x} [i, i] = z_{x i}

Z_{y} [i, i] = z_{y i}

W_{x} [i, i] = w_{x i}

W_{y} [i, i] = w_{y i}

Then its gradient is calculated as follows:

\frac{\partial {‖ u ‖}_{{DRTV}_{x}}}{\partial u} = C_{x}^{T} Z_{x} W_{x} C_{x} v_{u},

(26)

\frac{\partial {‖ u ‖}_{{DRTV}_{y}}}{\partial u} = C_{y}^{T} Z_{y} W_{y} C_{y} v_{u} .

(27)

DRTV-ASD-POCS algorithm

The pseudo-code for the proposed DRTV algorithm is presented in Algorithm 1. Specifically, in Algorithm 1, Lines 5 and 6 implement the POCS component, which ensures data consistency and enforces non-negativity constraints. The data consistency is achieved through algebraic reconstruction technique (ART)³⁷ operations combined with forward projection. Lines 13 to 17 correspond to the ASD component, which minimizes the DRTV norm. The gradient computation of the DRTV norm in Lines 14 and 15 is derived based on Equations (12) and (13), respectively.

Algorithm 1.

Pseudocodes of the DRTV-ASD-POCS Algorithm

1. Initialization:

β = 1.0; β_{red} = 0.995; n g = 4; α = 0.2; α_{red} = 0.95; r_{red} = 0.95; r_{\max} = 0.95

u = 0

3. repeat main loop

u_{0} = u

5. for

i = 1, N_{d}

u = u + β A_{i}^{T} \frac{g_{i} - A_{i} \times u}{A_{i}^{T} \times A_{i}}

POCS (ART)

6. for

i = 1, N_{i}

do if

u_{i} < 0

then

u_{i} = 0

Enforce Positivity

u_{r e s} = u

\tilde{g} = A u

d d = ‖ \tilde{g} - {\tilde{g}}_{0} ‖_{2}

10.

d p = ‖ u - u_{0} ‖_{2}

11. if {first iteration} then

d t v g = α * d p

12.

u_{0} = u

13. for

i = 1, n g

do DRTV-ASD loop

14.

d u_{x} = \nabla_{u} ‖ u ‖_{D R T V_{x}}; d u_{y} = \nabla_{u} ‖ u ‖_{D R T V_{y}}

15.

d u_{x} = d u_{x} / ‖ d u_{x} ‖_{2}; d u_{y} = d u_{y} / ‖ d u_{y} ‖_{2}

16.

u = u - d t v g * d u_{x} - b * d t v g * d u_{y}

17. end for

18.

d g = ‖ u - u_{0} ‖_{2}

19.

if d g > r_{m a x} * d p

and

d d > ε

then

d t v g = d t v g * α_{r e d}

20.

β = β * β_{r e d}

21. until {stopping criteria}

22. return

u_{r e s}

Reconstruction parameters

The complete specification of the DRTV algorithm consists of model parameters and algorithm parameters. The model parameters determine the solution of the optimization model, including the system matrix $A$ , projection method, gradient information $d u_{x}$ and $d u_{y}$ . The algorithm parameters usually do not affect the solution of the optimization model, but may affect the convergence behavior, rate and path of the algorithm, including $β$ , $β_{red}$ , $n g$ , $α$ , $α_{red}$ , $r_{red}$ , $r_{\max}$ and the steepest-descent step-size $d t v g$ are inherent to the algorithm and can be determined based on the Algorithm 1, and the weight b are discussed in Section 4.

Results

Correctness and convergence analysis of algorithm

In this section, the correctness of the DRTV algorithm is verified by using the FORBILD phantom. The size of the phantom is 256 × 256, the rotation center is located at the position of [128,128] in the center of the image, the length of the detector bin element is 1, the number of detector bin is the same as the length of the phantom, and the projection data are uniformly collected in the range of [0, π] for 360 angles for reconstruction. The Siddon ray-driven method is used to find the system matrix A and generate projection data of size 256 × 360.Each pixel of the grayscale image can represent 256 different gray levels ranging from 0 to 255.When the difference in the grayscale of the image is small enough, the human eye cannot perceive the difference between the two images. In the case of sufficient and ideal experimental data, the difference between the reconstructed image and the real image can be measured by root mean square error (RMSE). In general, the monitor is no longer able to distinguish the reconstructed image from the true image, and thus the proposed model, solution algorithm, and computer implementation are correct when $R M S E (u_{n}, u_{truth}) \leq 1 0^{- 4}$ . RMSE is defined as follows:

R M S E (u_{n}, u_{truth}) = \frac{{‖ u_{n} - u_{truth} ‖}_{2}}{\sqrt{length (u)}},

(28)

where

u_{n}

is the reconstructed image vector to be solved,

u_{truth}

is the truth image vector, and

length (u)

is the length of the image vector.

Figure 1 (a) shows the true image of the FORBILD phantom, and (b) shows the reconstructed image. It can be seen that (a) and (b) are almost identical and it is difficult to distinguish them with the naked eye. (c) shows a comparison of the vertical centerline profile between the true image and reconstructed image, and (d) shows a comparison of the horizontal centerline profiles between the true image and reconstructed image. It can be seen from the figure that the centerline profile almost completely overlaps. It is shown that the algorithm achieves a highly accurate reconstruction and meets the correctness metric for qualitative observations in this experiment.

Figure 1.

The algorithm correctness verification for FORBILD phantom (a) True image of FORBILD phantom; (b) Reconstructed image; (c) Comparison of vertical centerline profile between the true image and the reconstructed image: the location of the blue line in (a); (d) Comparison of horizontal centerline profile between the true image and the reconstructed image: the location of the red line in (a).

We introduce three metrics to evaluate the convergence behavior of the algorithm by observing the changes in these metrics. They are defined as follows:

M_{1} (n) = R M S E (u_{n}, u_{truth}),

(29)

M_{2} (n) = ‖ g - A u_{n} ‖_{2},

(30)

M_{3} (n) = \frac{| {‖ u_{n} ‖}_{T V} - {‖ u_{truth} ‖}_{T V} |}{{‖ u_{truth} ‖}_{T V}} .

(31)

$M_{1} (n)$ denotes the root-mean-square error between the reconstructed image and the true image; $M_{2} (n)$ denotes the $ℓ_{2}$ norm between the projection of the reconstructed image and the original projection data; $M_{3} (n)$ denotes the relative error between the TV value of the reconstructed image and the TV value of the true image. $u_{n}$ denotes the reconstructed image of the nth iteration and $u_{truth}$ denotes the true image.

Figure 2 (a) shows the iterative trend of RMSE for FORBILD phantom. Above 100 iterations, their RMSEs are less than 10⁻⁴, which meets the metric for quantitatively analyzing the experimental correctness. In summary, the qualitative observation and quantitative analysis of the reconstruction results can prove that the verification of the correctness of the algorithm is successful. (b) shows the iterative trend of the data error, and it can be seen that the trend is still decreasing. (c) shows the iterative trend of the relative error of their TV values, respectively. Since the TV value is the sum of the gradient magnitude variations of all pixel points of the image, assigning the TV relative error to each pixel point, the relative error of each pixel point is less than 10⁻⁴, the monitor has been unable to distinguish the changes of the image. But there is an upward dithering during the reconstruction process. This phenomenon is a characteristic of the TV algorithms in the iterative process, which often appears in the iterative process. It does not affect the convergence of this algorithm. After the vibration, its relative error will continue to decrease, which indicates that the solution of the optimization model still tends to converge. Therefore, the three image quality metrics all have reached the allowed state of convergence.

Figure 2.

Convergence analysis of reconstruction results for the FORBILD phantom (a) $M_{1} (n)$ iterative trend; (b) $M_{2} (n)$ iterative trend; (c) $M_{3} (n)$ iterative trend.

Evaluation of sparse reconstruction ability

In order to evaluate the sparse reconstruction ability of the DRTV algorithm, in this section, the experiments are performed at 20, 30, 40, and 50 projection angles for FORBILD phantom and real CT image by using TV, DTV, RTV and DRTV algorithm. The sizes of FORBILD phantom and real CT image are 256 × 256. In the experimental analysis, RMSE and structural similarity (SSIM) are used as the metrics of image reconstruction quality.

S S I M (x, y) = \frac{(2 μ_{x} μ_{y} + c_{1}) (2 σ_{x y} + c_{2})}{(μ_{x}^{2} + μ_{y}^{2} + c_{1}) (σ_{x}^{2} + σ_{y}^{2} + c_{2})},

(32)

where x is the truth image, y is the reconstructed image;

μ_{x}

and

μ_{y}

are the mean values of x, y respectively;

σ_{x}^{2}

and

σ_{y}^{2}

are the variances of x, y respectively;

σ_{x y}

is the covariance of x and y;

c_{1}

and

c_{2}

are very small constants. RMSE denotes the root-mean-square error of the pixel values of the two images, which reflects the error level of the two images. In addition, the more similar the two images are, the closer the value of SSIM is to 1.

Figure 3 (a) and (b) shows the comparison of the RMSE of the FORBILD phantom and the real CT image by using TV, DTV, RTV and DRTV algorithm for 20 projection angles, respectively. The DRTV algorithm achieves superior reconstruction accuracy after 500 iterations. Figure 4 shows the reconstruction results of the four algorithms for FORBILD phantom. At 20 projection angles, the reconstruction results of the TV and DTV algorithm have obvious artifacts and noise, and several image details are also smoothed out. From 20 to 50 projection angles, the accuracy of the reconstructed image increases. However, the reconstructed image still has some obvious artifacts compared with the original image at 50 projection angles. The reconstructed images by using the RTV and DRTV algorithms are almost indistinguishable from the original images with the naked eye at 20 projection angles. As the projection angle increases, the reconstruction results are still consistent with the original image. Figure 5 shows the reconstructed images on the real CT image by using the four algorithms. Compared to the DRTV algorithm, the TV and RTV algorithm clearly smooth out more image details, while the DTV algorithm has artifacts and noise, especially in the experimental results at 20 angles. The DRTV algorithm exhibits superior capability in preserving structural integrity and fine details of the reconstructed images.

Figure 3.

The comparison of RMSE at 20 sparse projection angles by using TV, DTV, RTV and DRTV algorithm for different phantom (a) FORBILD phantom; (b) real CT image.

Figure 4.

The comparison of the TV, DTV, RTV and DRTV algorithm for reconstruction of the FORBILD phantom: the number above the image indicates the number of projections; the text on the left indicates the algorithm used.

Figure 5.

The comparison of TV, DTV, RTV and DRTV algorithm for reconstruction of real CT image: the number above the image indicates the number of projections; the text on the left indicates the algorithm used.

Table 1 and Table 2 show the RMSE and SSIM of the reconstruction of the FORBILD phantom and real CT image using the TV, DTV, RTV and DRTV algorithm for 20 to 50 projection angles, respectively. DRTV algorithm's final convergence accuracy is 0.0116 and 0.0849, respectively. However, the final convergence accuracies of TV, DTV and RTV algorithm are only 0.0504, 0.0187, 0.0179 and 0.0903, 0.0880, 0.0890, respectively. It is clear that the DRTV algorithm provides higher reconstruction accuracy and structural similarity for the reconstruction.

Table 1.

The comparison of RMSE and SSIM for the FORBILD phantom by the TV, DTV, RTV and DRTV algorithm.

Number of projections		20	30	40	50
RMSE	TV	0.0504	0.0264	0.0122	0.0045
	DTV	0.0187	0.0150	0.0073	0.0011
	RTV	0.0179	2.6624e-04	1.7108e-04	1.2624e-04
	DRTV	0.0116	4.2656e-05	2.7713e-05	1.8905e-05
SSIM	TV	0.9877	0.9963	0.9988	0.9994
	DTV	0.9979	0.9985	0.9992	0.9995
	RTV	0.9980	0.9995	0.9995	0.9995
	DRTV	0.9981	0.9995	0.9995	0.9995

Table 2.

The comparison of RMSE and SSIM for real CT image by the TV, DTV, RTV and DRTV algorithm.

Number of projections		20	30	40	50
RMSE	TV	0.0903	0.0681	0.0594	0.0518
	DTV	0.0880	0.0669	0.0570	0.0513
	RTV	0.0890	0.0659	0.0578	0.0511
	DRTV	0.0849	0.0644	0.0558	0.0501
SSIM	TV	0.9498	0.9716	0.9784	0.9836
	DTV	0.9520	0.9723	0.9801	0.9840
	RTV	0.9520	0.9739	0.9798	0.9842
	DRTV	0.9561	0.9749	0.9812	0.9847

Stability analysis of algorithm

To evaluate the stability of the DRTV algorithm, the experiments are performed at 50 projection angles with different levels of Gaussian noise added to each projection data by using the TV, DTV, RTV and DRTV algorithm for FORBILD phantom and real CT image, respectively. The mean value of Gaussian noise is 0, and the variance is 0.01, 0.02, 0.03, 0.04 and 0.05, respectively. In the experimental analysis, RMSE and SSIM are used as the metrics of image reconstruction quality.

We select the reconstruction results of the FORBILD phantom and real CT image at 50 projection angles and add Gaussian white noise with variance of 0.05 to the projection data for elaboration. As shown in Figure 6 (a) to (d) and Figure 7 (a) to (d)are the reconstruction images of TV, DTV, RTV and DRTV algorithm, respectively, and (e) to (h) are the local area magnification images of the respective reconstruction results. Figure 8 shows the comparison of the RMSE by using TV, DTV, RTV and DRTV algorithm for the FORBILD phantom and the real CT image at 50 projection angles, and with the addition of Gaussian white noise with variance of 0.05 to the projection data, respectively. It can be seen that the TV, DTV, RTV and DRTV algorithm all show noise suppression properties. However, the comparison reveals that the reconstruction results of the DRTV algorithm are more stable and clearer. In Table 3, the reconstruction accuracies of these four algorithms on the FORBILD phantom are 0.0145, 0.0122, 0.0111, and 0.0044, respectively. In Table 4, the reconstruction accuracies of these four algorithms on the real CT image are 0.0764, 0.0657, 0.0636, and 0.0573, respectively. The quantitative analysis shows that the DRTV algorithm is superior in terms of stability and protection of image structure.

Figure 6.

(a)∼(d) The reconstruction of the FORBILD phantom results using TV, DTV, RTV, and DRTV algorithm under the condition of adding Gaussian white noise with variance of 0.05 to the projection data at 50 projection angles; (e)∼(h) Local area magnification images of the reconstruction results.

Figure 7.

(a)∼(d) The reconstruction of real CT image results using TV, DTV, RTV, and DRTV algorithm under the condition of adding Gaussian white noise with variance of 0.05 to the projection data at 50 projection angles; (e)∼(h) Local area magnification images of the reconstruction results.

Figure 8.

The comparison of the RMSE using TV, DTV, RTV and DRTV algorithm for different phantom at 50 projection angles and adding Gaussian white noise with a variance of 0.05 to the projection data (a) FORBILD phantom; (b) real CT image.

Table 3.

The comparison of RMSE and SSIM for FORBILD phantom under different levels of noise by TV, DTV, RTV and DRTV algorithm.

Noise level		0.01	0.02	0.03	0.04	0.05
RMSE	TV	0.0080	0.0102	0.0116	0.0129	0.0145
	DTV	0.0066	0.0085	0.0100	0.0113	0.0122
	RTV	0.0097	0.0106	0.0100	0.0108	0.0111
	DRTV	0.0027	0.0028	0.0039	0.0043	0.0044
SSIM	TV	0.9992	0.9990	0.9989	0.9987	0.9985
	DTV	0.9993	0.9992	0.9990	0.9989	0.9988
	RTV	0.9991	0.9990	0.9990	0.9989	0.9989
	DRTV	0.9994	0.9994	0.9994	0.9994	0.9994

Table 4.

The comparison of RMSE and SSIM for real CT image under different levels of noise by TV, DTV, RTV and DRTV algorithm.

Noise level		0.01	0.02	0.03	0.04	0.05
RMSE	TV	0.0746	0.0773	0.0751	0.0746	0.0764
	DTV	0.0650	0.0652	0.0653	0.0654	0.0657
	RTV	0.0631	0.0629	0.0632	0.0637	0.0636
	DRTV	0.0566	0.0571	0.0570	0.0569	0.0573
SSIM	TV	0.9654	0.9627	0.9649	0.9654	0.9637
	DTV	0.9743	0.9741	0.9740	0.9739	0.9737
	RTV	0.9760	0.9761	0.9759	0.9755	0.9755
	DRTV	0.9806	0.9803	0.9804	0.9804	0.9802

Discussions

The selection of the parameter b has a significant effect on the reconstruction results. As an adjustment factor, b is used to balance the weights in the x and y directions during the reconstruction process. A suitable value of b can ensure that the details of the image in different directions are reasonably preserved, thus improving the reconstruction quality. If b is not properly selected, it may lead to the loss of information or artifacts in a certain direction, thus affecting the overall reconstruction effect. When the value of b is small, which prefers the x direction, the reconstruction process may focus more on keeping the details of the image in the horizontal direction, resulting in sharper horizontal edges. However, if the value of b is too small, it may trigger excessive preservation in the x direction, introducing noise or artifacts. When the value of b is big, which prefers the y direction. For images with distinct directional features, the selection of b value can help to highlight features in a particular direction. Adjusting b can target the enhancement of structures in a certain direction, thus improving the performance of the reconstruction results. Therefore, optimizing the selection of b is crucial for improving reconstruction accuracy and stability.

To better understand the effect of parameter b, we performed additional experiments across a continuous range of b values. Figure 9 illustrates the reconstruction of a real CT image at 50 projection angles for different values of b (0.7, 0.8, 0.9, 1.0, and 1.1). It can be observed that the reconstruction quality improves with b = 0.9, which retains the most image details and provides sharper contours. The values of RMSE and SSIM in Table 5 further confirm that the image reconstruction quality is optimal when b = 0.9, with minimal error and high similarity to the original image. As b moves away from 0.9, either towards smaller or larger values, the reconstruction quality starts to degrade.

Figure 9.

(a)∼(e) The comparison of different b value (0.7, 0.8, 0.9, 1.0, 1.1) for reconstruction of real CT image. (f)∼(j) Local area magnification images of the reconstruction results.

Table 5.

The comparison of RMSE and SSIM for real CT image by with different b value.

b	0.7	0.8	0.9	1.0	1.1
RMSE	0.0535	0.0534	0.0501	0.0511	0.0546
SSIM	0.9827	0.9827	0.9847	0.9842	0.9819

In summary, optimizing b is crucial for enhancing reconstruction accuracy and stability. The optimal value of b is context-dependent and varies based on the structure of the object and the noise characteristics in the measurement data. For isotropic objects, b≈1 is typically sufficient, while for anisotropic objects, tuning b to emphasize the dominant direction can significantly improve reconstruction results. We recommend that the value of b be carefully selected based on the specific features of the object to achieve the best possible reconstruction outcome.

Conclusions

In this paper, a directional relative TV algorithm for sparse-view CT reconstruction is proposed, which combines the advantages of RTV with those of DTV, designed under the ASD-POCS framework. The algorithm effectively protects the structural features and texture details of the image during the reconstruction process. Experiments conducted on simulated phantoms and real CT images demonstrate the correctness, convergence and superior performance of the DRTV algorithm in sparse-view reconstruction. Compared with the TV, DTV, and RTV algorithm, the DRTV algorithm exhibits superior preservation of structural features and texture details. The qualitative and quantitative results show that the algorithm can realize high-precision image reconstruction while suppressing artifacts. Moreover, the approach is applicable to other medical imaging modalities.

Footnotes

Acknowledgements

This work was supported in part by National Natural Science Foundation of China under grant 62071281, and by Local Science and Technology Development Fund Project Guided by the Central Government under grant YDZJSX2021A003.

ORCID iDs

Yanan Wang

Yu Wang

Peng Liu

Chenyun Fang

Yanjun Zhang

Ruotong Yang

Zhiwei Qiao

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Pan

Sidky

Vannier

. Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction? Inverse Prob 2009; 25: 123009.

Zou

Sidky

, et al. Region of interest reconstruction from truncated data in circular cone-beam CT. IEEE Trans Med Imaging 2006; 25: 869–881.

Donoho

. Compressed sensing. IEEE Trans Inf Theory 2006; 52: 1289–1306.

Candes

Romberg

Tao

. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory 2006; 52: 489–509.

Sidky

Kao

Pan

, et al. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT. J X-ray Sci Technol 2006; 14: 119–139.

Sidky

Pan

. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol 2008; 53: 4777–4807.

Karl

Kalra

, et al. Clinical low dose CT image reconstruction using high-order total variation techniques. Phys Med Imaging 2010; 7622: 76225D.

Qiao

Wang

, et al. Study of CT image reconstruction algorithm based on high order total variation. Optik (Stuttg) 2020; 204: 163814.

Zhou

, et al. Adaptive-weighted high order TV algorithm for sparse-view CT reconstruction. Med Phys 2023; 50: 5568–5584.

10.

Tian

Jia

Yuan

, et al. Low-dose CT reconstruction via edge-preserving total variation regularization. Phys Med Biol 2011; 56: 5949–5967.

11.

Liu

Fan

, et al. Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction. Phys Med Biol 2012; 57: 7923–7956.

12.

Yan

Xia

, et al. Structure extraction from texture via relative total variation. ACM Trans Graph 2012; 31: 1–10.

13.

Chen

Jin

, et al. A limited-angle CT reconstruction method based on anisotropic TV minimization. Phys Med Biol 2013; 58: 2119–2141.

14.

Ning

Wei

. An algorithm for image reconstruction based on lp norm. Acta Phys Sin 2013; 62: 174212.

15.

Chambolle

Pock

. An introduction to continuous optimization for imaging. Acta Numer 2016; 25: 161–319.

16.

Chambolle

Pock

. A first-order primal-dual algorithm for convex problems with applications to imaging. J Math Imaging Vis 2011; 40: 120–145.

17.

Sidky

Jørgensen

Pan

. Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm. Phys Med Biol 2012; 57: 3065–3091.

18.

Sidky

Chartrand

Boone

, et al. Constrained TpV minimization for enhanced exploitation of gradient sparsity: application to CT image reconstruction. IEEE J Transl Eng Health Med 2014; 2: 1800418.

19.

Rigie

Riviere

. Joint reconstruction of multi-channel, spectral CT data via constrained total nuclear variation minimization. Phys Med Biol 2015; 60: 1741–1762.

20.

Wang

Nakamoto

Zhang

, et al. Reweighted anisotropic total variation minimization for limited-angle CT reconstruction. IEEE Trans Nucl Sci 2017; 64: 2742–2760.

21.

Zhang

Chen

Xia

, et al. Directional-TV algorithm for image reconstruction from limited-angular-range data. Med Image Anal 2021; 70: 102030.

22.

Qiao

. A simple and fast ASD-POCS algorithm for image reconstruction. J Xray Sci Technol 2021; 29: 491–506.

23.

Epel

Kotecha

Halpern

. In vivo preclinical cancer and tissue engineering applications of absolute oxygen imaging using pulse EPR. J Magn Reson 2017; 280: 149–157.

24.

Qiao

Zhang

Pan

, et al. Optimization-based image reconstruction from sparsely sampled data in electron paramagnetic resonance imaging. J Magn Reson 2018; 294: 24–34.

25.

Qiao

Liang

Tang

, et al. Optimization-Based image reconstruction from fast-scanned, noisy projections in EPR imaging. IEEE Access 2019; 7: 19590–19601.

26.

Qiao

Redler

Epel

, et al. A balanced total-variation-Chambolle-Pock algorithm for EPR imaging. J Magn Reson 2021; 328: 107009.

27.

Qiao

. An iterative reconstruction algorithm based on total nuclear variation for multi-channel EPRI. Optik (Stuttg) 2023; 287: 171114.

28.

Liu

Fang

Qiao

. An edge-preserving total nuclear variation minimization algorithm in EPR image reconstruction. Biomed Signal Process Control 2024; 87: 105426.

29.

Qiao

Liu

Fang

, et al. Directional TV algorithm for image reconstruction from sparse-view projections in EPR imaging. Phys Med Biol 2024; 69: 115051.

30.

Fang

Epel

, et al. Directional TV algorithm for fast EPR imaging. J Magn Reson 2024; 361: 107652.

31.

Liu

Xiong

Yang

, et al. A generalized relative total variation method for image smoothing. Multimed Tools Appl 2016; 75: 7909–7930.

32.

Zeng

Chen

, et al. Bilateral weighted relative total variation for low-dose CT reconstruction. J Digit Imaging 2023; 36: 458–467.

33.

Gong

Zeng

. Adaptive iterative reconstruction based on relative total variation for low-intensity computed tomography. Signal Processing 2019; 165: 149–162.

34.

Fei

Wei

Xiao

. Iterative directional total variation refinement for compressive sensing image reconstruction. IEEE Signal Process Lett 2013; 20: 1070–1073.

35.

Zhao

Pan

, et al. Sparse-view CT reconstruction based on gradient directional total variation. Meas Sci Technol 2019; 30: 055404.

36.

Siddon

. Fast calculation of the exact radiological path for a three-dimensional CT array. Med Phys 1985; 12: 252–255.

37.

Gordon

Bender

Herman

. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J Theor Biol 1970; 29: 471–481.