Variable Weighted Ordered Subset Image Reconstruction Algorithm

We propose two variable weighted iterative reconstruction algorithms (VW-ART and VW-OS-SART) to improve the algebraic reconstruction technique (ART) and simultaneous algebraic reconstruction technique (SART) and establish their convergence. In the two algorithms, the weighting varies with the geometrical direction of the ray. Experimental results with both numerical simulation and real CT data demonstrate that the VW-ART has a significant improvement in the quality of reconstructed images over ART and OS-SART. Moreover, both VW-ART and VW-OS-SART are more promising in convergence speed than the ART and SART, respectively.


INTRODUCTION
Many image reconstruction problems can be modeled by the following linear system: where y = (y 1 , y 2 , . . . , y n ) T ∈ R n is the observed data, n is the number of the projection datum, x = (x 1 , x 2 , . . . , x m ) T ∈ R m is the original image (T denotes the transpose of a matrix, R is the real number field, and m is the number of pixel, namely, image array is √ m × √ m) and R = (r i j ) is an n × m matrix. Image reconstruction is to estimate the original image x from the observed data y. Since (1) is usually ill-conditioned and the data in practice are noisy and very large, a solution by direct method involving the matrix R is infeasible [1]. Instead, some iterative methods, such as the algebraic reconstruction technique (ART) [2], simultaneous algebraic reconstruction technique (SART) [3], and expectation maximization (EM) algorithms [4], have been developed for image reconstruction. The iterative algorithm is a promising approach to achieve a better image quality in CT. However, limitations exist with respect to the required computation time.
In order to accelerate the convergence of iterative algorithms in image reconstruction, several improvements were designed based on the SART and EM algorithms. One remarkable technique is the ordered subset (OS) (also called block iterative (BI)) versions of the simultaneous schemes [5]. There are some guidelines to be considered [6] in the selection of subsets and blocks, but the relaxation parameters are fixed in numerical implementations [7,8]. In this paper, we propose variable weighted OS-SART (VW-OS-SART) and variable weighted ART (VW-ART) algorithms. Since the weights are considered in these subsets, the reconstructed images by VW-ART are better than those by OS-SART, VW-OS-SART, and ART. The reconstruction speeds of these algorithms are faster than those by OS-SART and those of ART.
The structure of the manuscript is as follows. In Section 2, we introduce a method for subset partition and corresponding weights selection. In Section 3, we derive the VW-OS-SART and VW-ART algorithms. In Section 4, we discuss the convergence property of the proposed algorithms. In Section 5, we apply our algorithms to simulated and practical data. In Section 6, we discuss some relevant issues and conclude the paper.

SUBSETS PARTITION AND WEIGHTS SELECTION
The main idea of ordered subset method is to partition the dataset into several subsets. It is crucial to find an effective way for the partitioning. In each iteration, we hope that the image is improved as much as possible with a lower computational cost. Hence, each subset should contain as 2 International Journal of Biomedical Imaging much complementary tomographic information as possible and have a relatively small size for computational balance. One of the natural subset partitions is formulated by dividing the projection data according to the projection angles. For example, for system (1), the 2D image array is There are 2n/T √ m the observed projection angles in each subinterval.
The projections are sequenced according to geometrical direction of the ray. The projection data of rays are sequenced according to geometrical position. Denote the index of projection datum by B = {1, 2, . . . , n}. Then we partition the index set B into T nonempty subsets Each B t contains all the indices of projection data in the tth angle subinterval in (2). The subsets B t are not necessarily disjoint.
The weights of the all projection data in one projection angle are equal. We define where α is the projection angle of y i . See Figure 1 for the geometrical interpretation of the weights selection. In each projecting angle subinterval, the weights locate on the two-segment line. Therefore the projection data index set B is divided into T subsets {B 1 , B 2 , . . . , B T }, and the weight corresponding to each subset is given by μ Bt (y i ) which satisfies Figure 1: Geometrical interpretation of the weights selection in formula (4), α denote projection angle, μ t (α) denote the weighted of the projection data at α projection angle in tth subinterval,

VARIABLE WEIGHTED VERSION OF OS-SART AND VARIABLE WEIGHTED ART ALGORITHMS
Let V and W be two positive definite diagonal matrices of order m and n, respectively. The following general Landweber scheme was studied in [1]: where λ l > 0 is the relaxation coefficient in the lth iteration. When then (6) is the SART [5] x (l+1) where R i is the ith row of R.
According to the partitioning and weighting strategies proposed in Section 2, the variable weighted version of the OS-SART algorithm (VW-OS-SART) can be written as where [l] = l mod (T) + 1 for l ≥ 0. The iteration process from l = νT to l = (ν + 1)T is called one cycle, in which each subset B t is visited exactly once. When then (6) is the same as (9). (6) becomes the ART [1]. The variable weighted version of this algorithm, called VW-ART, can be written as where i = l mod (n) + 1, and {μ Bt (y i )} (2ν + 1)n ≤ l < 2(ν + 1)n, ν = 0, 1, 2, . . . .
where R i and W i are the ith row of R and the ith diagonal element of W, respectively. Denote (2) . . . where m dimension vector, then we have a linear algebraic system Equation (1) is consistent if and only if (15) is consistent, and they have exactly the same solution set. The OS version of the Landweber scheme of (15) can be written as where (16) is the same as the VW-OS-SART formula (9).
Similar to that in [5], we get the following convergence results. Theorem 1. If 0 < λ l < 2, for all l ≥ 0, l λ l = ∞, and lim l→∞ λ l = 0, then the sequence {x (l) } generated by (9) converges to x ( * ) + P V x (0) even with inconsistent data. Here, P V is the orthogonal projection from R m to the kernel of R with respect to the inner product ·, · V .
VW-ART algorithm is a special ART algorithm, where the relaxation parameters {λ l } are special sequence, namely, the relaxation parameters are λ l μ i (where μ i ≤ 1). If the sequence {x (l) } generated by the ART algorithm converges, the sequence{x (l) } generated by the VW-ART algorithm (11) converges, too. That is the following theorem. Theorem 2. If 0 < λ l μ l < 2, for all l ≥ 0, l λ l = ∞, and lim l→∞ λ l = 0, then the sequence {x (l) } generated by (11) converges to x ( * ) + P V x (0) .
In order to study the speed of convergence, we call finishing one VW cycle 2 iteration numbers, because the consuming time of one VW-OS-SART cycle is almost twice as that of one OS-SART cycle (T/2 ordered subsets), and the time of one VW-ART cycle is twice as that of one ART cycle [9]. Namely, the iteration numbers in variable weighted algorithm are always even.

EXPERIMENT RESULTS
To test the performance of the VW-OS-SART, and compare it with OS-SART and VW-ART, we code all of them in the programming language C on a PC (RAM: 256 MB, CPU: Intel P4 2.8 GHz) under Windows XP operating system. The image array is 256 × 256, with 256 gray levels. Projections are simulated (or obtained in experiment) over 360 • with the uniform sampling scheme. There are totally 360 projections, 256 identical detector bins per projection. Pseudo random noises satisfying Gauss distribution are added to the Shepp-Logan phantom data, with variance being 2.5% of image mean, and the relaxation sequence λ l = 2. The numbers of subsets T = 72, each B t containing nine projections, that is, m = 256 2 , n = 256 × 360. The weights are computed as in Section 1. For example, if y is one of the 256 projections with the projection degree 1 • , then μ A1 (y) = 0.2, μ A72 (y) = 0.8; if y is one of the 256 projections with the projection degree 5 • , then μ B1 (y) = 1; if y is a projection with the projection degree 8 • , then μ B1 (y) = 0.4, μ B37 (y) = 0.6; if y is one of the 256 projections with the projection degree 10 • , then μ B37 (y) = 1.
In Figure 2, we present a Shepp-Logan phantom and the images reconstructed by the above algorithms (the initial value of iterative reconstruction is x (0) = (0, 0, . . . , 0)). It can be seen from Figure 2 that all the reconstructed images by VW-ART are better than those by OS-SART, VW-OS-SART (relaxation parameter is 2) and ART. The reconstructed speeds of VW-OS-SART are faster than those by OS-SART; the reconstruction speeds of VW-ART are faster than those of ART (see Figure 3).
Since the true image is known in the Shepp-Logan phantom experiment, to quantify the performance of the VW-OS-SART and VW-ART formulas, the mean square error is computed from the images in the CT phantom experiment. In Figure 3, we present the difference of the mean square error of the reconstructed images by ART, VW-ART, OS-SART, and VW-OS-SART.
In the practical data experiments we use a cylindrical object of 80 mm in diameter, with 14 apertures of different diameters inside, and scan with X-ray generator with 220 KVp and 10 mA. The whole scan time is 1.5 second. The distance between the focus and the detector is 1.000 mm; the distance between the center of rotation and the detector is 50 mm. Figure 4  The images reconstructed by the VW-ART (Figure 4(e)) are better than those reconstructed by the OS-SART, the VW-OS-SART, and the ART (Figures 4(b)-4(d)) in quality. The criterion of comparison is the resolution of images, namely, whether the apertures can be detected.

DISCUSSIONS AND CONCLUSION
In our experiments, we have found that for both VW-OS-SART and VW-ART, the partition of ordered subset affects the quality of reconstruction image. And the selection of weights for each projection is crucial to variable weighted algorithms, and it can be a piece-wise linear function, Gaussian function and so on. The convergence conditions for VW-OS-SART and VW-ART are guidelines for the selection of the relaxation parameters. The relaxation parameters in ART algorithm and VW-ART algorithm are generally smaller than those in SART algorithm and VW-OS-SART algorithm. The partition of the subsets is generally based on the geometrical relationship of projections. The weights selection is given by linearity in this paper, and it can be given by an exponential as well, which is similar to ordered subsets method, Jinxiao Pan et al. generally, "increasing the number of the subsets accelerates iterative convergence, but there is a point beyond which image quality degrades due to lack of either tomographic or statistical information within subsets" [6]. The number of subsets may affect the convergence of the algorithm. Instead of partitioning the data in the projection domain, another way is to partition the pixel set in the image domain, which is under investigation.
In conclusion, we have proposed the VW-OS-SART and VW-ART algorithms and established their convergence. The proposed algorithms have been evaluated with Shepp-Logan phantom and practical data. In VW-OS-SART and VW-ART algorithms, the sum of the weighted of each projection data is 1, the weighting varies with the geometrical direction of the ray, the weight and partitioning procedures are intuitively reasonable and easy to be implemented. We use a weighted partitioning scheme such that each subset contains equal-spaced angle X-rays only, then these algorithms can have a significant improvement in the quality of reconstructed images over ART and SART, the experimental results on both digital phantom and real CT data have demonstrated that the VW-ART has a significant improvement in the quality of reconstructed images over ART and SART; both VW-ART and VW-OS-SART are more promising in convergence speed than ART and SART, respectively.

PROOF OF THEOREM 1
Let x ( * ) be the minimal V norm solution of the problem, where the V norm is defined as x 2 V = x, x V = Vx, x , and P V is the orthogonal projection from R m to the kernel of R (the kernel of C is identical with that of R) with respect to the inner product ·, · V , and R V ,W is the operator norm from R m (with ·, · V ) to R n (with ·, · W ). The following lemma is proved in [1].
International Journal of Biomedical Imaging Because (1) is consistent if and only if (15) is consistent, the formula (16) is the same as formula (9). The following lemma of Theorem 1 is proved in [1]. then the sequence {x (l) } generated by (9) converges to x ( * ) + P V x (0) , even if (1) is inconsistent.
To establish the convergence of the proposed algorithms, we estimate C V ,W with the corresponding V and W as follows: where the inequality is according to the convexity of the function t → t 2 . Therefore, C V ,W ≤ 1, we can choose ρ to be 1, and get the convergence result.