Fast parallel implementation for total variation constrained algebraic reconstruction technique

2.1 TV constrained ART algorithm

In this study, we focus on the reconstruction of the two-dimensional case. The CT reconstruction problem can be formulated as a linear equation system Af = p(1) where f∈R ^N is the unknown image to be reconstructed, p∈R ^M is the measured projection data, and A = (a _ij ) is an M×N system matrix, whose element a _ij , named weight coefficient, represents the contribution of the jth pixel to the ith ray integral. In practice, the computation of a _ij is determined by the system model, which has critical effect on both the accuracy and the efficiency of the reconstruction. The most used system model is the line intersection model. In this model, the element a _ij is calculated as the intersection length of the jth pixel and the ith ray. The ART algorithm solves Eq. (1) by the following iterative procedures: f j ( w + 1 ) = f j ( w ) + λ p i - ∑ n = 1 N a in f n ( w ) ∑ n = 1 N a in 2 a ij , w = 0 , 1 , … i = w mod M + 1(2) where w is the iterative number, λ is a relaxation parameter, and f⁽⁰⁾ is an initial estimate of the solution. In most cases, f⁽⁰⁾ is initialized to a zero vector.

To improve the reconstruction quality of ART from sparse-view projections, the TV constraint is usually incorporated into its reconstruction procedure. Thus, the reconstruction algorithm of TV constrained ART can be formulated as follows: min f ∥ f ∥ TV , s . t . Af = p(3) where ∥f ∥ _TV represents the TV norm of image f. Here, the TV of the reconstructed image is defined as the L₁ norm of its gradient image, which can be calculated as: ∥ f ∥ TV = ∥ ∇ f ∥ 1 = ∑ s , t ( f s , t - f s - 1 , t ) 2 + ( f s , t - f s , t - 1 ) 2(4) where ∇ is the local gradient operator, s and t are, respectively, the row and column index of image f. To solve Eq. (3), the first step is to perform a whole iterative reconstruction using Eq. (2). Then, the gradient descent method is used to minimize the TV of image f, which is given by f ( w , m + 1 ) = f ( w , m ) − α d ( w ) v ‖ v ‖ 2 d ( w ) = ∑ s , t ( f s , t ( w , 0 ) − f s , t ( w − 1 , N T V − 1 ) ) 2(5) where m is the iterative number of TV minimization during the wth ART iteration, N _TV is the maximum number of iterations for each ART iteration, a is the relaxation factor, and v is the partial derivative of ∥f ∥ _TV. The element v _{s ,t} of v is calculated as: v s , t = ∂ ∥ f ∥ TV ∂ f s , t ≈ 2 f s , t - f s - 1 , t - f s , t - 1 ( f s , t - f s - 1 , t ) 2 + ( f s , t - f s , t - 1 ) 2 + ɛ - f s + 1 , t - f s , t ( f s + 1 , t - f s , t ) 2 + ( f s + 1 , t - f s + 1 , t - 1 ) 2 + ɛ - f s , t + 1 - f s , t ( f s , t + 1 - f s , t ) 2 + ( f s , t + 1 - f s - 1 , t + 1 ) 2 + ɛ(6)

where ɛ is a small positive number.

The above two steps are implemented alternately until the termination condition is satisfied. Obviously, to improve the reconstruction speed of ART-TV algorithm, we should try to accelerate both the ART iteration and the TV minimization. Since the system matrix plays a crucial role during the iterative procedure of ART, we will firstly propose a fast algorithm to calculate the system matrix. Then we utilize the parallel computing technique to further accelerate the two steps. To perform the forward-projection and back-projection operations, we use two arrays, a and b, to store the pixel indices and their corresponding intersection lengths. In addition, the total number of elements is recorded into a variable n₁.

The image f to be reconstructed can be represented by a square grid (i.e., reconstruction region), which consists of N = n×n pixels. Let f _j be the image value of the jth pixel, where j is the pixel index, 0≤j < N. The ray SE is considered as a line segment connecting the X-ray focal spot S(x _S , y _S ) and the center E(x _E , y _E ) of a detection bin, as shown in Fig. 1. The equation of the ray SE can be written as g(x, y) = Ax + By + C = 0, where A = y _E −y _S , B = x _S −x _E , and C = y _S ×x _E −x _S ×y _E .

OPEN IN VIEWER

Fig. 1

Reconstruction region and the arrangement of the image pixels.

First, we should determine whether the ray passes through the square grid. Obviously, when a ray passes through the reconstruction region, it will intersect with at least one of the two diagonals of the square grid, as shown by the dashed line in Fig. 1. For example, suppose SE intersects with diagonal PD, then the vertex P and D will lie on the two sides of SE respectively. In this case, we have g(−q, q)×g(q, −q) < 0. Likewise, if g(−q, −q)×g(q, q) < 0, then the ray intersects with the other diagonal.

When the ray SE passes through the reconstruction region, we need to determine the two intersection points of the ray with the sides of the square region and calculate their coordinates and pixel indices. Let L and R be the left and right intersection point, whose pixel indices are n _L and n _R , respectively. Let (x _L , y _L ) and (x _R , y _R ) denote their coordinates, here we suppose x _L ≤ x _R . If x _E  = x _S , then the ray is perpendicular to the horizontal axis, and it will pass through n pixels. The pixel index n _L within the first row can be calculated by n _L  =⌊(x _L  + q)/δ⌋. The pixel indices of other rows can be obtained by adding a constant integer n to n _L in turn. The intersection lengths of all these pixels are δ. Let k denote the slope of the ray SE. If x _E ≠ x _S , then k = (y _E − y _S ) / (x _E − x _S ). Let y _L denote the y-coordinate of the intersection point of SE with line x =−q, where y _L  = (A×q−C)/B. If |y _L | < q, then point L intersects with the left side of the reconstruction region, where x _L  =−q, y _L  = y _L , and n _L  = n×⌊(q− y _L )/δ⌋. If y _L  > q, then point L intersects with the top side of the reconstruction region, where x _L  =−(B×q + C)/A, y _L  = q, and n _L  =⌊(x _L  + q)/δ⌋. If y _L  < -q, then point L intersects with the bottom side of the reconstruction region, where x _L  = (B×q−C)/A, y _L  =−q, and n _L  =⌊(x _L  + q)/δ⌋+ n×(n−1). The coordinates (x _R , y _R ) of point R, as well as its corresponding pixel index n _R can be calculated in the same way as point L.

Without loss of generality, we consider the case of 0≤k≤1 in the following discussion. After calculating the intersection information of points L and R, we begin to calculate the pixel indices and the corresponding intersection lengths of those pixels that intersect with the ray SE. For ease of calculation, we denote the total length intercepted by a column as l _x , and that intercepted by a row as l _y , where l x = δ 1 + k 2 l y = δ 1 + k 2 / | k |(7)

Next, we consider that the ray SE passes through only one row. If ⌊y _L /δ⌋=⌊y _R /δ⌋, then SE intersects with only one row. Let u and l denote the index and intersection length of the first pixel in a row, respectively. If point L is on the bottom side of the reconstruction region, then the intersection length l of pixel n _L may be equal to l _x , or less than l _x , as shown in Fig. 2.

OPEN IN VIEWER

Fig. 2

The ray passes through only one row. (a) shows the intersection lengths of the pixels are equal to l _x . (b) shows the intersection length of the last pixel of the upper ray and that of the first pixel of the lower ray are less than l _x .

For 0≤k ≤1, we have l = l _x − (x _L /δ -⌊x _L /δ⌋)×l _x . Then, for other pixels intersecting with SE in the same row, their pixel indices and intersection lengths can be easily calculated. If point R is on the top side of the reconstruction region, we can refer to the calculation of point L. To determine when to stop our algorithm, we define a remainder length denoted by l _r . Once we have stored the intersection length of one or more pixels, the remainder length will be updated by subtracting the stored length. For this purpose, we calculate the length l₀ from point R to point L and initialize l _r to the value of l₀. The pseudocode of the ray SE (0≤k≤1) passing through only one row is presented in Algorithm 1.

If ⌊y _L /δ⌋ ≠ ⌊y _R /δ⌋, the ray will pass through several rows of the square grid. We first determine the principal direction of SE according to its slope k. If |k|,> 1, then the principal direction is along the x-axis, otherwise, it is along the y-axis. For the case of 0≤k ≤1, the pixels are processed row by row. Let N _y denote the maximum number of pixels whose intersection length is l _x within the same row, and r _y denote the remainder length except these N _y pixels, where N _y  =⌊l _y /l _x ⌋, r _y  = l _y − N _y ×l _x . When SE passes through a row, there are usually the following cases, as shown in Fig. 3.

OPEN IN VIEWER

Fig. 3

Different cases for the ray passing through a row. (a) shows the ray passes through N _y  + 1 pixels, where r _y  < l < l _x . (b) shows the ray passes through N _y  + 2 pixels, where 0 < l < r _y . (c) shows the ray passes through N _y  + 1 pixels, where l = l _x . (d) shows the ray passes through N _y  + 1 pixels, where l = r _y .

If r _y  < l < l _x , as Fig. 3(a) shows, then the ray SE will pass through N _y  + 1 pixels. The index of the last pixel in the row is u + N _y , whose intersection length is r _y  + l _x −l. The intersection length of those pixels between u and u + N _y are l _x . If 0 < l < r _y , as Fig. 3(b) shows, then the ray SE will pass through N _y  + 2 pixels. The index of the last pixel in the row is u + N _y  + 1, whose intersection length is r _y −l. The intersection length of those pixels between u and u + N _y  + 1 are l _x . If l = l _x or l = r _y , as Fig. 3(c) and 3(d) shows, then the ray SE will pass through N _y  + 1 pixels. If r _y ≠ 0, the ray SE will pass through N _y  + 1 pixels and intersect with either the lower left vertex of the first pixel, or the upper right vertex of the last pixel. The pixel index of the last pixel in the row is u + N _y . If r _y  = 0, then SE will pass through N _y pixels and intersect with both the lower left vertex of pixel u and the upper right vertex of pixel u + N _y −1. The intersection lengths of those pixels from u to u + N _y −1 are l _x .

According to the index and the intersection length of the last pixel, we can easily determine the intersection information of the first pixel of the next row to be processed. Hence, the above procedure will be repeated until l _r is equal to zero. The pseudocode of the proposed algorithm for a given ray SE (0≤k≤1) passing through several rows is presented in Algorithm 2.

As mentioned above, the performance of ART-TV algorithm includes two steps: ART iteration and TV minimization. For ART, the iterative procedures are implemented in a particular sequence of projection views. Therefore, it is not appropriate for parallel implementation in general. However, for a given projection view, the reconstruction could be implemented in parallel providing that none of the pixels is updated by two or more rays simultaneously. Thus, we can decompose the reconstruction task of a given projection view into several subtasks. For simplicity, we just decompose the reconstruction task evenly by the index of detector bins. To avoid the access collision, the iteration of each subtask is performed according to the ascending or descending order of its bin indices. Then, we use multithreading technique to accelerate these subtasks. Open multi-processing (OpenMP) is a parallel programming model for shared memory multi-core platform. The OpenMP application programming interface (API) specifies a set of compiler directives and a library of subroutines, which enable us to easily parallelize the sequential code [30,31, 30,31]. In this work, we use the OpenMP to parallelize the implementation of ART. Let T _n be the number of threads, and B _n be the number of detector bins. Then, the pseudocode of the multithreading-based ART is presented in Algorithm 3.

In recent years, GPU parallel computing technique has been widely used in various fields of scientific computing [32–35]. The compute unified device architecture (CUDA), a C-like API, provides an easy way to develop parallel program for GPU. As discussed above, the TV minimization consists of four steps: computing d^(w), computing the elements v _{s ,t} of the partial derivative v, normalizing v, and updating image f. Therefore, we design three CUDA kernels, kernel_1, kernel_2 and kernel_3, to perform the TV minimization. The kernel_1 is to compute the partial derivative image v using Equation (6), kernel_2 is to normalize the partial derivative image, i.e. v = v/∥ v ∥ ₂, and kernel_3 is to update image f using Equation (5). The pseudocode of the GPU accelerated TV minimization is presented in Algorithm 4.

3.1 Efficiency of our proposed algorithm

We simulated a fan-beam CT system with a linear detector of 2048 bins spaced by 0.192 mm, 1150 mm source to detector distance, 650 mm source to origin distance, and 180 projection views over 360°. The two-dimensional FORBILD head phantom was used in the numerical experiments, and the simulated projection data were generated using analytical method from this phantom. To evaluate the efficiency of our method, we reconstructed the phantom using ART-TV algorithm with image sizes of 1024×1024 and 2048×2048. Correspondingly, the pixel sizes were 0.209×0.209 mm² and 0.1045×0.1045 mm², respectively. For comparison, we also used the Siddon method and the distance-driven method to calculate the system matrices. Table 1 shows the computation time of the three methods with one iteration.

For 1024×1024 images, as shown in Table 1, the forward projection using our method is 4.8 times and 1.7 times faster than that using the Siddon and the distance-driven methods, respectively. While for 2048×2048 images, our method is 4.3 times and 2.0 times faster than the Siddon and the distance-driven methods, respectively. Overall, for the reconstruction of 2048×2048 images, our method can obtain speedup factors of 2.6x and 1.6x compared with the Siddon and the distance-driven method, respectively. The results show that our method is very efficient compared with the Siddon method. It should be noted that the back-projection operation is fully determined by the elements of the system matrix calculated during the forward projection operation. As can be seen from Table 1, the computation time of the back-projection using Siddon’s method is almost the same as our proposed method for a given image.

The Siddon method contains two main time-consuming processes, i.e. merging the intersection parametric sets into an ordered set and calculating the pixel indices, which take up about 26% and 41% of the total computational time, respectively. These processes involve multiple operations, such as multiplication, division and integer rounding, which decrease the efficiency of Siddon’s method greatly. However, in our method, most of the pixel indices are calculated using only addition or subtraction operations. Besides, most intersection lengths are calculated by assigning a constant value l _x or l _y for a given ray. Therefore, our method can save a considerable amount of computation time compared with the Siddon method.

To evaluate the accuracy of our method, we performed the ART-TV reconstruction using the Siddon method and our method, respectively. The reconstruction image was 2048×2048, whose values were initialized to zero. The relaxation factor λ was set to 0.2, and the projections were accessed in a random permutation scheme. For the total variation, the regularization parameter α was set to 0.26, ɛ was set to 10^–8, and the iteration number was 20. For comparison, we also reconstructed the phantom using the distance-driven method and the analytical FBP method. Note that, the Shepp-Logan filter and linear interpolation were used in the reconstruction of FBP. Figure 4 shows the reconstructed images using different methods.

OPEN IN VIEWER

Fig. 4

The original phantom and the reconstructed images. (a) shows the original phantom. (b) and (c) show the results using Siddon’s method after 10 and 80 iterations, respectively. (d) and (e) show the results using distance-driven method after 10 and 80 iterations, respectively. (f) and (g) show the results using our method after 10 and 80 iterations, respectively. (h) shows the result using FBP algorithm.

As can be seen from Fig. 4, the visual quality of the reconstructions using Siddon’s method, distance-driven method and our method is similar. With the increasing of the iteration number, the reconstruction quality is getting better. Because the projection data are incomplete for the analytical FBP algorithm, severe artifacts can be found in Fig. 4(h). To evaluate our method in a quantitative way, the normalized root mean square error (NRMS) and the normalized mean absolute error (NMA) were computed in the experiment. The formulas of the two errors are given by NRMS = ∑ j = 1 N [ f j - r j ] 2 ∑ j = 1 N [ f j - f ¯ ] 2(8) NMA = ∑ j = 1 N | f j - r j | ∑ j = 1 N | f j |(9) where f _j is the pixel value of the phantom, r _j is the reconstructed pixel value, f ¯ is the mean value of the phantom, N is the total number of pixels. For this purpose, we discretized the two-dimensional FORBILD head phantom into 2048×2048 pixels of equal size. Table 2 lists the values of the two errors of Fig. 4.

As can be seen from Table 2, the NRMS and NMA errors of our method are almost identical with those of Siddon’s method, which indicates that our method is very accurate in calculating the system matrix. Both errors of the distance-driven method are smaller than those of the LIM-based method for the same number of iterations, which demonstrates the superiority of the distance-driven method in iterative reconstruction. With the number of iterations increasing, both errors are becoming smaller. We can obtain satisfactory reconstruction results after large numbers of iterations. Moreover, the results show that the reconstruction quality of ART-TV method is much better than that of the conventional FBP algorithm for sparse-view projections.

In this experiment, different numbers of threads were used to accelerate the ART reconstruction for a 2048×2048 image. Meanwhile, the GPU-based method was used to accelerate the TV minimization, in which the size of the CUDA block was set to 16 ×16, and the grid was set to 128×128. Figure 5 shows the computation time of the parallel reconstruction after one iteration using the Siddon and our method.

OPEN IN VIEWER

Fig. 5

Computation time of the parallel reconstruction using the Siddon method and our method.

As can be seen from Fig. 5, the total reconstruction time using the Siddon method with one thread is 46.2 seconds, and that using the proposed method is 15.2 seconds. By using GPU implementation, the computation time of the TV minimization is reduced to less than 0.8 seconds. With the number of threads increasing, the computation time of the forward and backprojection process is gradually reduced. Using ten threads, the total reconstruction time of the proposed method is reduced to 2.2 seconds. We obtained a speedup factor of 23x compared with the conventional CPU implementation. To assess the reconstruction quality of the parallel reconstruction, Tables 3 and 4 show the NRMS and NMA errors of the two methods after 10 iterations with different numbers of threads.

From Tables 3 and 4, we can see that our parallel reconstruction has almost the same accuracy as the conventional single-thread CPU reconstruction. The results show that our parallel implementation is efficient and accurate. Figure 6 shows the reconstructed images with 10 threads after 10 iterations and their profiles comparison of the two methods.

OPEN IN VIEWER

Fig. 6

Reconstructed images and their profiles. (a) shows the parallel reconstruction using Siddon’s method. (b) shows the parallel reconstruction using our proposed method. (c) shows the central vertical profiles of the original phantom and the reconstructed images (a) and (b).