|
 
 |
| ORIGINAL ARTICLE |
|
| Year : 2020 | Volume
: 10
| Issue : 1 | Page : 1-11 |
|
Computationally efficient system matrix calculation techniques in computed tomography iterative reconstruction
Golshan Mahmoudi1, Mohammad Reza Ay1, Arman Rahmim2, Hossein Ghadiri1
1 Department of Medical Physics and Biomedical Engineering, Tehran University of Medical Sciences; Research Center for Molecular and Cellular Imaging, Tehran University of Medical Sciences, Tehran, Iran
2 Department of Radiology and Physics, University of British Columbia; Department of Integrative Oncology, BC Cancer Research Centre, Vancouver, BC, Canada
| Date of Submission | 10-Jun-2019 |
| Date of Decision | 27-Jul-2019 |
| Date of Acceptance | 04-Sep-2019 |
| Date of Web Publication | 06-Feb-2020 |
Correspondence Address:
Dr. Hossein Ghadiri
Department of Medical Physics and Biomedical Engineering, Tehran University of Medical Sciences, Tehran
Iran
 Source of Support: None, Conflict of Interest: None
DOI: 10.4103/jmss.JMSS_29_19
Background: Relative to classical methods in computed tomography, iterative reconstruction techniques enable significantly improved image qualities and/or lowered patient doses. However, the computational speed is a major concern for these iterative techniques. In the present study, we present a method for fast system matrix calculation based on the line integral model (LIM) to speed up the computations without compromising the image quality. In addition, we develop a hybrid line–area integral model (AIM) that highlights the advantages of both LIM and AIMs. Methods: The contributing detectors for a given pixel and a given projection view, and the length of corresponding intersection lines with pixels, are calculated using our proposed algorithm. For the hybrid method, the respective narrow-angle fan beam was modeled by multiple equally spaced lines. The computed system matrix was evaluated in the context of reconstruction using the simultaneous algebraic reconstruction technique (SART) as well as maximum likelihood expectation maximization (MLEM). Results: The proposed LIM offers a considerable reduction in calculation times compared to the standard Siddon algorithm: 2.9 times faster. Differences in root mean square error and peak signal-to-noise ratio were not significant between the proposed LIM and the Siddon algorithm for both SART and MLEM reconstruction methods (P > 0.05). Meanwhile, the proposed hybrid method resulted in significantly improved image qualities relative to LIM and the Siddon algorithm (P < 0.05), though computations were 4.9 times more intensive than the proposed LIM. Conclusion: We have proposed two fast algorithms to calculate the system matrix. The first is based on LIM and was faster than the Siddon algorithm, with matched image quality, whereas the second method is a hybrid LIM–AIM that achieves significantly improved images though with its computational requirements.
Keywords: Area integral model, computed tomography, forward and back projection, iterative image reconstruction, line integral model, system matrix
How to cite this article:
Mahmoudi G, Ay MR, Rahmim A, Ghadiri H. Computationally efficient system matrix calculation techniques in computed tomography iterative reconstruction. J Med Signals Sens 2020;10:1-11 |
How to cite this URL:
Mahmoudi G, Ay MR, Rahmim A, Ghadiri H. Computationally efficient system matrix calculation techniques in computed tomography iterative reconstruction. J Med Signals Sens [serial online] 2020 [cited 2023 Jun 5];10:1-11. Available from: https://www.jmssjournal.net/text.asp?2020/10/1/1/277822 |
| Introduction | |  |
Given the inherent risk of ionizing radiation, reduction of imaging doses in computed tomography (CT) imaging is highly desirable.[1],[2],[3],[4] Today, iterative reconstruction techniques are one of the most important strategies for reducing radiation dose in CT. Indeed, compared to conventional methods such as filtered back projection, iterative reconstruction techniques have significant potential to this end.[4],[5],[6],[7],[8] Despite their advantages, the critical disadvantage of iterative methods is the computational speed of image reconstruction.[9]
Let us formulate the CT image reconstruction problem in a linear algebra framework:
AX = P (1)
Where X is the unknown image of N = n2 pixels, Pis the measured projection data of M = v*d rays (v = number of projection views and d = number of detector elements), and A is the system matrix (M × N) that relates the image space to the projection space. The columns of the system matrix correspond to the pixels of the reconstructed image, and the rows are the measured projections. Hence, each element, aij (sometimes called weighting coefficient), of the matrix represents the contribution of pixel j to ray integral i. The system matrix has large dimensions, and working with this very large matrix challenges the reconstruction task. Moreover, an immense amount of memory may be required to save the system matrix.[10],[11] To tackle this, computations are commonly performed on the fly, involving forward and back projections (which are transpose operations of one another). These operations are the most time-consuming components in iterative reconstruction algorithms.[9],[11],[12],[13],[14],[15],[16]
Two different models have been used to perform the projection operations on CT systems: line integral model (LIM) and area integral model (AIM). In the former, the detector is considered as a point and so the rays are effective of zero width, connecting the source to the center of the detector cell. One way to do LIM is to compute the intersection length between the jth pixel with the ith ray. Another approach is to do bilinear or trilinear interpolation in the line model which is like calculating the perpendicular distance between the center of the pixel and the line.[17] For the LIM approach, an effective method was proposed by Siddon.[18] In this algorithm, a series of horizontal and vertical lines form a grid. Then, the intersection points of a line (ray) with grid limits are calculated. Two arrays of intersection points are formed: one for intersection with the horizontal lines and the other for the vertical lines. The mentioned two arrays are then merged into one array and sorted, which are used to calculate the intersection length and pixel indices.
In the AIM, the rays are considered as a narrow-angle fan beam that connects the source to the detector cell and computes the intersection area between the jth pixel and the ith ray.[19] The LIM has low computational complexity but lower image quality due to undersampling and aliasing problems. The AIM is much closer to a real CT system and overcomes the sampling problems but has high computational complexity.
In the present study, we present a method to speed up the computation of projection operations of the system matrix without compromising the image quality. Here, we considered the two-dimensional (2D) fan-beam CT geometry and developed our method based on LIM. However, our work could be extended to finite-size beam and 3D cone-beam CT geometry. Furthermore, using this method, we developed a hybrid LIM–AIM that approximates the area integral by averaging multiple line integrals. Finally, to evaluate the performance of the new methods in terms of both speed and image qualities, the computed system matrix was used to reconstruct images of a numerical phantom, using the simultaneous algebraic reconstruction technique (SART) and the maximum likelihood expectation maximization (MLEM) for both noiseless and noisy data. We are going to implement the proposed hybrid method in the context of micro-CT in the next study.
| Methods | | |
System matrix calculation method
For simplicity of presentation, we assumed a nonoverlapping uniform pixel model and fan-beam geometry with zero width line beam that connects the X-ray (point) source to the center of the detector cell. The idea is to consider only the contributing detectors at each projection view. More elaborately put, for a given pixel and a given projection view, the contribution from the pixel only comes from the detectors that are bounded by the beams passing through the entire pixel [Figure 1]. To implement this in the algorithm, we need to calculate two angles. The first angle α is the angle between the lines that pass through pixel boundaries such that cover the whole pixel. The second angle β is the angle between the line pass through the pixel center and the line which connects the source to the first detector cell [Angles α and β have been shown in [Figure 1]. To get this angle, we need to know the pixel center coordinates. As focal spot rotates around the object, the pixel center coordinates change with rotation. As such, it needs to calculate the new pixel center coordinates in the new coordinate system at each projection view. In the original coordinate system that is considered here, the detector is perpendicular to the y-axis. Finally, the new pixel center coordinates are: | Figure 1: Schematic representation of the proposed algorithm: (a) All projection rays for a pixel and a given angle of projection, (b) the contributing detectors are bounded by the beams passing through the entire pixel
Click here to view |
Where xrot and yrot are the new pixel center coordinates and xcent and ycent are the pixel center coordinates in the original coordinate system and θ is the projection view. Now, we can find the lower and upper index bounds of the contributing detectors. The lower index bound is determined by subtracting half of the angle α from the angle β and the upper index bound by adding instead of subtracting. Therefore, the lower and upper index bounds of the contributing detectors are:
Where arc_det is the angle of the detector cell's arc and is calculated by dividing the fan angle to the number of detector elements. After finding the contributing detectors, which are often only a few, the length of the intersection line from these detectors with the pixel should be calculated. At last, we can calculate the weighting coefficients (system matrix elements) by normalizing the intersection lengths.
The above calculation was related to the curve detectors, as shown in [Figure 1]. We can develop our method for the flat-panel detectors. The lower and upper index bounds of the contributing detectors for the flat-panel detectors are:
Where size_pix is the pixel size, size_det is the detector element size, and d is the number of detector elements. SID and SDD are the source to isocenter distance and source to detector distance, respectively.
As mentioned above, for the Siddon algorithm, the intersection points of a ray with horizontal and vertical lines (the pixel grid) are calculated separately and form two intersection arrays. This is needed to merge these two arrays and sort them to calculate the intersection length. In addition, multiple operations such as multiplication, division, and integer rounding should be performed to calculate the pixel indices. On the contrary, in our method, the intersection length between a pixel and a ray was calculated for a given pixel and a given projection view, so no intersection point arrays are formed to require merging and sorting. Furthermore, in the presented algorithm, merely addition and subtraction operations are used to calculate the pixel indices which are simpler than operations used in the Siddon algorithm. Therefore, the presented algorithm is faster than the Siddon algorithm mainly due to its simplicity and decreased computational effort. In addition, this is a general method to speed up system matrix calculation by considering the contributing detectors. As such, unlike Siddon algorithm, it could be easily generalized to the fan-beam geometry with a finite-beam size and could be applied to other methods based on AIM. Moreover, the proposed algorithm could be extended to the 3D geometry and is highly parallelizable.
As the major part of the system matrix consists of zeroes, it is possible to store this matrix as a sparse matrix. In a sparse matrix, only the nonzero elements and its indices are stored. Thus, significantly less storage is required, and the data processing will proceed faster. For this reason, we define the system matrix as a sparse matrix in the algorithm (equally applicable to our method as well as Siddon method).
Hybrid line–area integral model
In the AIM, the intersection area of a narrow-angle fan beam with pixels is calculated, which is an improved approximation, and results in better image quality but has more computational complexity. In the LIM, the intersection length of a beam with pixels is calculated, which is simpler but has lower accuracy. Therefore, we have developed a hybrid method that used the advantages of both LIM and AIM and approximates the area integral by averaging multiple line integrals. Line integral calculations to approximate the area integral have less computational complexity than the accurate calculation of area integral. On the other hand, approximating the area integral is more closely related to real CT imaging than the line integral. As such, that method has the potential to combine the advantages of LIM and AIM. More elaborately put, a narrow fan beam can be modeled by multiple equally spaced lines that connect the source to one detector cell [Figure 2]. In this way, the area integral can be approximated by averaging the line integrals considered for a narrow fan beam. The system matrix element can be obtained as follows: | Figure 2: Schematic representation of modeling the finite-size beam with multiple equally spaced lines in the proposed hybrid line–area integral model: (a) zero width beam, and (b) finite-size beam
Click here to view |
Where n is the number of lines per detector and li, j is the weighting coefficient of each line. The number of lines needed to model the beam is critical to the accuracy of the system matrix and therefore to image quality. Using more lines per detector is more closely related to reality. However, increasing the number of lines per detector increases the computation time. Furthermore, the number of lines per detector depends on the detector element size: for the large detector element, more lines need to be considered. Therefore, to choose the number of lines per detector, we provide a trade-off between computation time and image quality.
Image reconstruction method
To evaluate the performance of the proposed models, SART[20] and MLEM[21] were used to reconstruct the images. SART can be expressed as:
Where k is the iteration number and λ is the relaxation parameter which is critical to the quality and speed of the reconstruction (here, we consider λ to be 0.1).
MLEM as a statistical iterative method which could be modeled the noise can be obtained as:
Where k is the iteration number.
Image quality assessment
We used the 2D Shepp–Logan head phantom[22] [Figure 3] to create the simulated projection data. The projection data were generated from the phantom analytically, using CT simulator data.[23] The noisy transmission data can be simulated as follows. Given the mathematical Shepp–Logan phantom, the line integral  was computed along the projection path or ray i. By the Lambert–Beer law,  , and the knowledge of  in routine clinical studies, the mean Ii was calculated. Then, the Poisson noise was added to the mean Ii. After the noisy transmission datum  was simulated from the noise-free transmission datum, the corresponding noisy projection datum pi was obtained by the logarithm transform of the noisy transmission datum,  . For noisy data, five different Poisson noise levels were used such that from level 1 to level 5 noise increases. The computations are performed for a typical clinical CT scanner with fan-beam geometry. The source-to-isocenter distance was 540 mm, and source-to-detector distance was 950 mm. The projection view number was 720. For each projection view, 512 detector cells were distributed with equal distance, which defines a field of view of 250 mm. The detector element size was 1.8 mm. The pixel size of the reconstructed image was varied among 1.9, 0.98, and 0.49 mm; correspondingly, the reconstructed image was 128 × 128, 256 × 256, and 512 × 512, respectively. All the algorithms were implemented with MATLAB/R2015a and tested on a PC platform (Intel i7-6700 Processor, 8.0 GB memory, 3.4 GHz CPU). We measure the performance of the new algorithm in terms of time and reconstruction quality.
To evaluate the reconstructed image quality quantitatively, we used two of the most common metrics used to measure accuracy: root mean square error (RMSE) that detects large errors in a small number of pixels[24] and peak signal-to-noise ratio (PSNR).[25] The equations of these two metrics are as follows:
Where Ii,j and Ri,j are the pixel values of the original and the reconstructed images, respectively, and MAXI is the maximum possible pixel value of the image.
To evaluate the spatial resolution of the reconstructed images, the modulation transfer function (MTF) was obtained. To compute the MTF, we first obtained the edge-spread function (ESF) along the profile as indicated by the white line in [Figure 3]. Then, the line-spread function (LSF) was calculated by the derivation of ESF. Finally, the MTF was obtained by applying the Fourier transform to the LSF and normalized the MTF such that MTF (0) = 1.
| Results | | |
System matrix calculation method
To assess computational efficiency, we evaluated our proposed algorithm and the popular Siddon algorithm for the fan-beam geometry with different pixel sizes. [Figure 4] shows the computational performance of the two methods for different pixel sizes. From [Figure 4], we can see that for 128 × 128, 256 × 256, and 512 × 512 images, our method is, respectively, 3.5, 2.9, and 2.8 times as fast as the Siddon algorithm. | Figure 4: The computation time of different sizes of the system matrix using the line integral model algorithm versus Siddon algorithm
Click here to view |
[Table 1] presents quantitative comparisons of reconstructed images by two error metrics of RMSE and PSNR, using two reconstruction methods of SART and MLEM for both noiseless and noisy data. As we can see from the table, the proposed LIM produced similar error measures (both RMSE and PSNR) to Siddon. In other words, our proposed method obtained faster computations of the system matrix, with no compromise to image quality for two cases of noiseless and noisy data. | Table 1: Quantitative comparisons of the reconstructed images based on the Siddon algorithm, proposed line integral model, and hybrid methods
Click here to view |
Hybrid line–area integral model
For the hybrid algorithm, [Figure 5] shows the difference between the original image and the reconstructed image (using RMSE) versus the number of lines per detector. As we expected, modeling the beam with multiple equally spaced lines reduced the RMSE, so using hybrid method could improve image quality compared to the LIM. The figure also shows that the difference does not change dramatically using >5 lines per detector for a detector element size of 1.8 mm (i.e., almost 3 lines per 1 mm of the detector). As such, in our work, we chose 5 lines per detector to perform the hybrid method. | Figure 5: Root mean square error versus number of lines per detector. Modeling a narrow fan beam using the proposed hybrid line–area integral model, decreasing the root mean square error, and improving the image quality
Click here to view |
The results with our proposed hybrid algorithm are also presented in [Table 1], indicating further improved quantitative performance. As we can see from the table, the differences of RMSE and PSNR between the proposed hybrid model and two other methods (the Siddon algorithm and the proposed LIM) are more evident in SART compared to MLEM.
In [Figure 6], RMSE of the images reconstructed by SART (after 20 iterations for noiseless data and after 6 iterations for noisy data) is shown as a function of the number of iterations for three system matrix calculation methods and both noiseless and noisy data. [Figure 7] shows these results for the images reconstructed by MLEM (up to 50 iterations). As we can see in [Figure 6]b with the increase of iteration number, RMSE decreases and then increases due to the instability of the iterative reconstruction methods. We can also see that the hybrid method resulted in better image quality than two other methods during all iterations for both reconstruction methods. Moreover, we can see that with increase of the iteration number, the difference of the RMSE between three methods becomes more evident, and this is more obvious for SART compared to MLEM. The proposed LIM and the Siddon algorithm produced similar RMSE methods during all iterations for both reconstruction methods. | Figure 6: Root mean square error of the images reconstructed by simultaneous algebraic reconstruction technique versus the number of iteration, using Siddon algorithm, the proposed line integral model, and hybrid model: (a) Noiseless data up to 20 iterations and (b) noisy data up to 6 iterations
Click here to view |
 | Figure 7: Root mean square error of the images reconstructed by maximum likelihood expectation maximization versus the number of iterations (up to 50 iterations), using the Siddon algorithm, the proposed line integral model, and hybrid model: (a) Noiseless data and (b) noisy data
Click here to view |
The graphical representation of [Table 1] (for noisy data) is shown in [Figure 8] and [Figure 9]. [Figure 8] shows the means (bars) and standard deviations (error bars) of RMSE values of the images reconstructed by SART (up to 6 iterations) and MLEM (up to 50 iterations) at five different Poisson noise levels for three system matrix calculation methods. Using the pairwise Wilcoxon test, significant differences in RMSE were seen in the proposed hybrid model compared to the proposed LIM and the Siddon algorithm for both reconstruction methods (P < 0.05), whereas differences were not significant between the proposed LIM and the Siddon algorithm for both reconstruction methods (P > 0.05). The same results were obtained for PSNR values shown in [Figure 9]. | Figure 8: Root mean square error of the images reconstructed by simultaneous algebraic reconstruction technique (after 6 iterations) and maximum likelihood expectation maximization (after 50 iterations), using three system matrix calculation methods of Siddon algorithm, the proposed line integral model, and hybrid model. The bars show the mean value of root mean square error for five different noise levels, while the error bars indicate the standard deviation
Click here to view |
 | Figure 9: Peak signal-to-noise ratio of the images reconstructed by simultaneous algebraic reconstruction technique (after 6 iterations) and maximum likelihood expectation maximization (after 50 iterations), using three system matrix calculation methods of Siddon algorithm, the proposed line integral model, and hybrid model. The bars show the mean value of root mean square error for five different noise levels, while the error bars indicate the standard deviation
Click here to view |
To evaluate the spatial resolution of the reconstructed images, the MTF curves were utilized. The measured MTF curves of images reconstructed by SART (after 20 iterations for noiseless data and after 6 iterations for noisy data) are shown in [Figure 10] for three system matrix calculation methods and both noiseless and noisy data. [Figure 11] shows these results for the images reconstructed by MLEM (up to 50 iterations). As illustrated in this figure, the proposed hybrid model can yield a better spatial resolution than the LIM methods for both noiseless and noisy data and so could preserve more edge details. Moreover, the results demonstrate that the spatial resolution is almost the same for the proposed LIM and the Siddon algorithm. As we expected, the spatial resolution is better for noiseless data than the noisy data for both reconstruction algorithms (SART and MLEM), and this is more evident for SART that is more sensitive to the noise. | Figure 10: MTF curves from images reconstructed by simultaneous algebraic reconstruction technique, using Siddon algorithm, the proposed line integral model, and hybrid model: (a) noiseless data, after 20 iterations and (b) noisy data, after 6 iterations
Click here to view |
 | Figure 11: Modulation transfer function curves from images reconstructed by maximum likelihood expectation maximization (after 50 iterations), using Siddon algorithm, the proposed line integral model, and hybrid model: (a) Noiseless data and (b) noisy data
Click here to view |
[Figure 12] and [Figure 13] present visual comparisons of reconstructed images for the Shepp–Logan phantom using SART (up to 6 iterations) and MLEM (up to 50 iterations), respectively. The proposed LIM and the Siddon algorithm depict fluctuations in the uniform region. On the contrary, the proposed hybrid method is much closer to that of the original phantom. The reconstruction results of the three methods appear visually comparable, though quantitatively, as shown in [Table 1], the performance was best for our proposed hybrid method. | Figure 12: Visual comparisons of reconstructed images of the Shepp–Logan phantom using simultaneous algebraic reconstruction technique (after 6 iterations) for both noiseless data and five different noise levels (left to right), using: Siddon ( first row), the proposed line integral model (second row), and the proposed hybrid model (third row)
Click here to view |
 | Figure 13: Visual comparisons of reconstructed images of the Shepp–Logan phantom using maximum likelihood expectation maximization (after 50 iterations) for both noiseless data and five different noise levels (left to right), using: Siddon ( first row), the proposed line integral model (second row), and the proposed hybrid model (third row)
Click here to view |
The profiles along horizontal sharp edge as shown in [Figure 14] and [Figure 15] demonstrated that the reconstructed images using the proposed hybrid method have a better spatial resolution than that of the proposed LIM and the Siddon algorithm.
| Discussion | | |
In this work, two fast algorithms for computation of the system matrix are proposed: (a) a fast algorithm based on the line integral and the intersection length calculation which improves the computation time significantly relative to the popular Siddon algorithm, while not compromising quantitative accuracy and (b) a hybrid LIM–AIM that combines the advantages of both LIM and AIM for better image quality and less computation time.
The proposed LIM algorithm was based on considering only the contributing detectors at each projection view, such that for a given pixel and a given projection view, the contribution from the pixel only comes from the detectors that are bounded by the beams passing through the pixel. As such, we obtain the lower and upper index bounds of the contributing detectors. In a practical example for a typical clinical CT scan, we can consider just ten detectors instead of 1000 detectors at each projection view, and therefore, the computational speed is greatly increased. In this method, the system matrix was calculated pixel by pixel: for a given pixel, the contributions of the pixel for all projection rays were calculated, and so on for other pixels. In this way, the system matrix is completed column by column during the implementation of the algorithm.
For the Siddon algorithm, the intersection points of a ray with the pixel grid are calculated separately and form two intersection arrays which need to be merged and sorted to calculate the intersection length. As such, sorting and merging of arrays and the calculation of pixel indices involve multiple operations of multiplication, division, and integer rounding function that are relatively time-consuming. For our method, by contrast, the intersection length between a pixel and a ray was calculated for a given pixel and a given projection view, and therefore, no intersection point arrays are formed to require merging and sorting. Furthermore, in the presented algorithm, mere addition and subtraction operations are used to calculate the pixel indices which are simpler than operations used in the Siddon algorithm. The computational complexity for the Siddon algorithm for each beam is O (n) and the overall complexity is roughly O (n·v·d) (where n2 = the number of pixels, v = number of projection views, and d = number of detector elements). For our algorithm, by contrast, the computational complexity for each pixel is O (v) and the total complexity is O (n. n. v), showing our algorithm to have less computational complexity than the Siddon algorithm. Moreover, in terms of image quality, the quantitative and qualitative analyses showed that differences in RMSE were not significant between the proposed LIM and the Siddon algorithm. As a result, our implementation scheme for computing the system matrix made our method more computationally efficient than the Siddon algorithm while maintaining accuracy. Moreover, the proposed algorithm could be extended to finite-beam size and 3D geometries and is highly parallelizable.
We described and studied the hybrid LIM–AIM to construct a more realistic and accurate CT system model by simulating a narrow fan beam with a set of equally spaced lines to approximate the intersection areas by averaging the line integrals. As such, in this method, instead of calculating the area integral, multiple line integrals are calculated and averaged. Therefore, this algorithm does not have the difficulty of computing the area integral and requires less computational effort compared to AIM. In addition, approximating the intersection area instead of calculating the intersection length is more reflective of the real CT geometry/physics and could improve the image quality and image resolution compared to LIM. As increasing the number of lines per detector increased both image quality and computation time, a trade-off between computation time and image quality should be provided to choose the optimum number of lines per detector. Furthermore, the number of lines per detector depends on the detector element size. Modeling the finite-size beam with five equally spaced lines per detector decreases RMSE compared to one line per detector. However, the image quality does not improve significantly using more than five lines per detector in the case of 1.8 mm detector pixel size. This should be pointed out that the optimum number of lines per detector could change based on the detector pixel size and the acquisition system geometry. As such, significant improvements are obtained by our model, though increasing computation by a factor of nearly five compared to the LIM.
The system matrix can be precomputed, stored, and reused, using a sparse matrix, so that the time of the entire reconstruction process is reduced. Alternatively, it can be generated on the fly during the reconstruction process, thus minimizing the memory usage of the system.
Adding multiple Poisson noise levels to the sinogram, the proposed hybrid model still maintains its superiority to two other methods and even becomes more obvious. Therefore, for low-dose CT applications, it is more practical to use the proposed hybrid model instead of LIM-based methods. Moreover, the effect of using the proposed hybrid method on image quality was seen more in algebraic reconstruction techniques (ART) compared to statistical iterative reconstruction (SIR). As such, at low-level noises, the proposed hybrid-based ART might be preferable to the proposed hybrid-based SIR.
Due to instability of iterative methods, RMSE decreases and then increases with increasing iteration numbers, and for higher noise levels, this increment happens at lower iteration numbers. Higher error bars for SART relative to MLEM imply that instability of iterative methods is more evident for nonstatistical iterative methods compared to statistical methods due to noise modeling in statistical iterative reconstruction. However, MLEM results in oversmoothing and also the uncertainty of edge detection compared to SART. Therefore, this method is not appropriate for CT, especially when the high-resolution images are desired.
The potential of the proposed hybrid method to improve image quality becomes more evident at high iteration numbers. Specifically, with increasing iteration numbers, RMSE between the hybrid method and two other methods increases. At low iteration numbers, the low-frequency components are reconstructed, and at higher iteration numbers, the higher frequency components are reconstructed. The high-frequency components are representations of image details and therefore impact spatial resolution. As such, our proposed hybrid method is more effective for reconstruction of high-frequency content and improves the spatial resolution of the image compared to two other methods. The MTF curves confirm this result and indicated a better spatial resolution for the proposed hybrid model. The LIM has lower image quality due to undersampling and aliasing problems. On the contrary, the proposed hybrid model increases the MTF cutoff and therefore could overcome sampling problems and aliasing artifacts.
In recent studies,[26],[27] ray modeling was shown to give very noticeable improvements. Here, our results do not show noticeable visible improvements. More complex phantoms and studies might show larger differences, which remain to be studied.
| Conclusion | | |
We have presented two fast algorithms to calculate the system matrix in CT image reconstruction: a fast and accurate LIM and a hybrid LIM–AIM. Results have shown that the proposed LIM method is faster than the Siddon algorithm but maintains accuracy. The improved computational efficiency was mainly due to the simplicity in algorithmic operations. The hybrid method was shown to yield better image qualities than our LIM-based method. Furthermore, as this method could save a large amount of computation time and memory compared to AIM-based algorithms, it could be a good alternative for high-accuracy imaging.
Financial support and sponsorship
This study was financially supported by a grant from Tehran University of Medical Sciences (grant number: 33474).
Conflicts of interest
There are no conflicts of interest.
| Biographies | | |
Golshan Mahmoudi is currently PhD candidate in Medical Physics at Tehran University of Medical Sciences, Tehran, Iran. Her research interests include medical imaging and medical image reconstruction.
Email: [email protected]
Mohammad Reza Ay is Professor of Medical Physics at Tehran University of Medical Sciences, Faculty of Medicine, Department of Medical Physics and Biomedical Engineering. His main research field is molecular and cellular imaging and authored significant number of papers in this area.
Email: [email protected]
Arman Rahmim is Associate Professor of Radiology and Physics & Astronomy, Faculties of Medicine and Science at the University of British Columbia. His researchfield is quantitative imaging and molecular imaging and authored significant number of papers in this field.
Email: [email protected]
Hossein Ghadiri is Assistant Professor of Medical Physics at Tehran University of Medical Sciences, Faculty of Medicine, Department of Medical Physics and Biomedical Engineering. His main research field is advance x-ray imaging systems and authored significant number of papers in this field.
Email: [email protected]
| References | | |
| 1. | Rampinelli C, Origgi D, Vecchi V, Funicelli L, Raimondi S, Deak P, et al. Ultra-low-dose CT with model-based iterative reconstruction (MBIR): Detection of ground-glass nodules in an anthropomorphic phantom study. Radiol Med 2015;120:611-7.  |
| 2. | Brenner DJ, Hall EJ. Computed tomography – An increasing source of radiation exposure. N Engl J Med 2007;357:2277-84. |
| 3. | Sodickson A, Baeyens PF, Andriole KP, Prevedello LM, Nawfel RD, Hanson R, et al. Recurrent CT, cumulative radiation exposure, and associated radiation-induced cancer risks from CT of adults. Radiology 2009;251:175-84. |
| 4. | Pontana F, Henry S, Duhamel A, Faivre JB, Tacelli N, Pagniez J, et al. Impact of iterative reconstruction on the diagnosis of acute pulmonary embolism (PE) on reduced-dose chest CT angiograms. Eur Radiol 2015;25:1182-9. |
| 5. | McCollough CH, Primak AN, Braun N, Kofler J, Yu L, Christner J. Strategies for reducing radiation dose in CT. Radiol Clin North Am 2009;47:27-40. |
| 6. | Fessler JA, Sonka M, Fitzpatrick JM, editors. Statistical image reconstruction methods for transmission tomography. In: Handbook of Medical Imaging. Medical Image Processing and Analysis. Vol. 2. Bellingham, WA: SPIE; 2000. |
| 7. | Thibault JB, Sauer KD, Bouman CA, Hsieh J. A three-dimensional statistical approach to improved image quality for multislice helical CT. Med Phys 2007;34:4526-44. |
| 8. | Wang G, Yu H, De Man B. An outlook on x-ray CT research and development. Med Phys 2008;35:1051-64. |
| 9. | Long Y, Fessler JA, Balter JM. 3D forward and back-projection for X-ray CT using separable footprints. IEEE Trans Med Imaging 2010;29:1839-50. |
| 10. | Van Hemelryck Tessa WS, Maggie G, Joost BK, Jan S. The Implementation of Iterative Reconstruction Algorithms in MATLAB: Master's thesis, Department of Industrial Sciences and Technology, University College of Antwerp. Belgium; 2007. |
| 11. | Ha S, Li H, Mueller K, editors. Efficient area-based ray integration using summed area tables and regression models. Proceedings of the 4 th International Meeting on Image Formation in X-ray CT; 2016. |
| 12. | Arcadu F, Nilchian M, Studer A, Stampanoni M, Marone F. A forward regridding method with minimal oversampling for accurate and efficient iterative tomographic algorithms. IEEE Trans Image Process 2016;25:1207-18. |
| 13. | De Man B, Basu S, editors. Distance-Driven Projection and Backprojection. IEEE Nucl Sci Symp Conf Rec. IEEE; 2002. |
| 14. | Miao C, Liu B, Xu Q, Yu H. An improved distance-driven method for projection and backprojection. J Xray Sci Technol 2014;22:1-8. |
| 15. | Gao H. Fast parallel algorithms for the x-ray transform and its adjoint. Med Phys 2012;39:7110-20. |
| 16. | Zhang S, Zhang D, Gong H, Ghasemalizadeh O, Wang G, Cao G. Fast and accurate computation of system matrix for area integral model-based algebraic reconstruction technique. Opt Eng 2014;53:113101. |
| 17. | Rahmim A, Cheng JC, Blinder S, Camborde ML, Sossi V. Statistical dynamic image reconstruction in state-of-the-art high-resolution PET. Phys Med Biol 2005;50:4887-912. |
| 18. | Siddon RL. Fast calculation of the exact radiological path for a three-dimensional CT array. Med Phys 1985;12:252-5. |
| 19. | Yu H, Wang G. Finite detector based projection model for high spatial resolution. J Xray Sci Technol 2012;20:229-38. |
| 20. | Jiang M, Wang G. Convergence of the simultaneous algebraic reconstruction technique (SART). IEEE Trans Image Process 2003;12:957-61. |
| 21. | Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomogr 1984;8:306-16. |
| 22. | Shepp LA, Logan BF. The Fourier reconstruction of a head section. IEEE Trans Nucl Sci 1974;21:21-43. |
| 23. | Ghadiri H, Rahmim A, Shiran MB, Soltanian-Zadeh H, Ay MR, editors. A Fast and Hardware Mimicking Analytic CT Simulator. 2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (2013 NSS/MIC); 2013. |
| 24. | Herman GT. Fundamentals of Computerized Tomography: Image Reconstruction from Projections. New York: Springer Science and Business Media; 2009. |
| 25. | Huynh-Thu Q, Ghanbari M. Scope of validity of PSNR in image/video quality assessment. Electron Lett 2008;44:800-1. |
| 26. | Hofmann C, Knaup M, Kachelrieß M. Effects of ray profile modeling on resolution recovery in clinical CT. Med Phys 2014;41:021907. |
| 27. | De Man B, Basu S. Distance-driven projection and backprojection in three dimensions. Phys Med Biol 2004;49:2463-75. |
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9], [Figure 10], [Figure 11], [Figure 12], [Figure 13], [Figure 14], [Figure 15]
[Table 1]
|
Search Pubmed for