An iterative image reconstruction algorithm for SPECT

An iterative image reconstruction algorithm for SPECT 1 corresp first-author Corresponding author, jwzhao@nju.edu.cn Corresponding author, jwzhao@nju.edu.cn jwzhao@nju.edu.cn 1 Physics School, Nanjing University, Nanjing210093, China Physics School, Nanjing University, Nanjing210093, China Properties of two algorithms for iterative reconstruction of SPECT images, LS-MLEM and LS-OSEM,are studied and compared with the ML-EM algorithm in this paper. By using projection data of heavy-noise, their effectiveness in improving SPECT image quality is evaluated. A phantom with hot and cold lesion is used in the investigation. The reconstructed images using LS-MLEM or LS-OSEM show that there is not a rapid increase in image noise,and the "best" estimate is assuming that the reconstructed images satisfy the statistical model. The major advantage of using LS-MLEM or LS-OSEM algorithm in SPECT imaging is in their ability to accurately control for heavy-noise. And LS-OSEM algorithm obviously improves the convergence rate. Least squares Maximum likelihood Image reconstruction Higher-noise I. INTRODUCTION 1 In single photon emission computed tomography (SPECT) applications, the projection data of radiotracer activity are heavily contaminated by statistical noise. The noise-control methods of iterative image reconstruction and incorporating statistical models can provide better noise properties than the conventional filtered back-projection method. The maximum-likelihood expectation-maximization(ML-EM) algorithm and ordered-subset variant (OSEM) algorithm [ 1 1 , 2 2 ] are widely applied. However, the impact of statistical noise is a key problem in the image reconstructions, and it is crucial to estimate uncertainties of the reconstructed images [ 3 3 ]. In practical applications, firstly, the iterative algorithm is terminated before convergence to control the noise in reconstructed images, for higher-noise (low-count) specially. Only a few iterative algorithms with explicit multiplicative update equations can be applied. So many iterative algorithms, such as those based on principles of maximum-likelihood, penalized-likelihood and Bayesian maximum a posteriori (MAP) estimation, are applied to image reconstruction. However, these iterative algorithms are nonlinear, and a theoretical formulation for noise propagation in image reconstruction has been intractable. Whichever noise-control algorithm is used, at least one arbitrary parameter (often several) shall be chosen, but how to choose it (them) is a key issue. Secondly, those iterative algorithms are derived by assuming a low-noise (high-count) approximation, which means that the noise in the reconstructed images is much smaller than the mean images. Therefore, they become less accurate at higher-noise. As iteration number increases, noise variance in a reconstructed image tends to increase. So, a noise-controlling iterative algorithm in ML-EM reconstruction stops before the ML point is reached [ 4 4 , 5 5 ], and the final image obtained depends on the stopping point. For statistical estimation problems, the best estimate can also be determined using the least squares (LS) and weighted least squares (WLS) criterion. The iterative algorithms proposed for LS or WLS image reconstruction include the optimal step size (OSS) algorithms, the gradient descent (GD) algorithms, the conjugate gradient (CG) algorithms, etc. The primary difference among them is the way to determine the step direction. CG algorithms are used the most widely in SPECT image reconstruction [ 6 6 - 8 8 6-8 ]. The step direction depends on the previous step direction and the gradient directions. The advantage of CG algorithms is that each new step dose not spoils the work of previous steps. However, this is not the case for LS/WLS-CG algorithms with a pre-conditioner, i.e. a transformation matrix. Performance of the LS/WLS-CG algorithms depends heavily on the pre-conditioner chosen. In general, the convergence rate of the WLS-CG algorithm is about ten times that of the ML-EM algorithm. Also, at large iteration numbers, the WLS-CG exhibits a faster increase in image noise than the ML-EM algorithm. The WLS-CG algorithm uses a Gaussian statistical model instead of the more correct Poisson model used in ML-EM algorithms. So, the WLS-CG algorithms are less accurate at higher-noise (low-count), and they are not easy to incorporate the non-negativity constraint in the reconstruction, unlike ML-EM algorithms. According to standard least squares, the mathematical form of the LS criterion is an equation. Because of the equation elements contaminated by statistical noise and the large dimension of the equation, the close-form solut ions of the equation are not often used in image reconstruction directly. The statistical noise is a Poisson model at higher-noise. Therefore, it is proposed that the equation based on standard least squares should be estimated using the ML-EM algorithms. Termed as LS-MLEM algorithm in this paper, it is advantageous in avoiding negativity and a faster increase in image noise at large iteration numbers. Mainly, the LS-MLEM based on least squares and ML-EM satisfies the Poisson model and a low-noise approximation. II. LS-MLEM ALGORITHM 1 In SPECT, the process of obtaining projection data is as follows: <math display="block" id="M1"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>C</mi></mstyle><mi>λ</mi><mn>,</mn></mrow></math> where P is the vector with the elements representing photon counts measured in a given projection bin, λ is the vector with the elements representing photon counts emitted from a given voxel, and C is a system matrix with the elements c il representing the relative contribution of l th voxel to projection bin p i . The matrix C may model the effects of attenuation, detector response, and/or scatter, but these effects are ignored in this paper, so as to evaluate the result of reconstructed image based on LS-MLEM algorithm. Since the measurement of P is contaminated by noise, P is frequently modeled as a random process, and the elements are independent Poisson random variables. The mathematical form of the LS criterion is given by <math display="block" id="M2"><mrow><mi>E</mi><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>=</mo><mo>|</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>C</mi></mstyle><mi>λ</mi><msup><mrow><mo>|</mo><mo>|</mo></mrow><mn>2</mn></msup><mn>.</mn></mrow></math> The squared residual error E ( λ ) shall be minimized. According to standard least squares, we have <math display="block" id="M3"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>C</mi></mstyle><mtext>T</mtext></msup><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo>=</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>C</mi></mstyle><mtext>T</mtext></msup><mstyle mathvariant="bold" mathsize="normal"><mi>C</mi></mstyle><mi>λ</mi><mn>,</mn></mrow></math> where C T is transpose matrix of the system matrix C . Defining X = C T P , its elements x j are determined by <math display="block" id="M4"><mrow><msub><mi>X</mi><mi>j</mi></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><msub><mi>c</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>p</mi><mi>i</mi></msub><mtext> </mtext><mtext> </mtext><mi>j</mi><mo>=</mo><mn>1,2</mn><mn>,</mn><mo>⋯</mo><mn>,</mn><mi>J</mi><mn>,</mn></mrow></math> where, I is the total number of given projection bin, and J is the total number of given voxel in a slice. P includes a Poisson noise. Therefore x j in Eq. (4) is a linear superposition of the projection data with Poisson noise. Therefore, the vector X , with its elements x j , also include Poisson noise. The close-form solution of Eq. (3) does not satisfy non-negativity constraint in the reconstruction. Signal-to-noise ratio SNR x of the vector X is given by SNR x ∝ I 1/2 SNR p , SNR p is signal-to-noise ratio of the vector P . Since the total number of projection bin is several thousands, the SNR x being ten times more than the SNR p can be approached. In Eq. (3), let H = C T C be a constant matrix based on system matrix, its elements h jl is determined by <math display="block" id="M5"><mrow><msub><mi>h</mi><mrow><mi>j</mi><mi>l</mi></mrow></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><msub><mi>c</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>c</mi><mrow><mi>i</mi><mi>l</mi></mrow></msub><mtext> </mtext><mtext> </mtext><mi>j</mi><mn>,</mn><mi>l</mi><mo>=</mo><mn>1,2</mn><mn>,</mn><mo>⋯</mo><mn>,</mn><mi>J</mi><mn>.</mn></mrow></math> Because Eq. (3) is of lower-noise, ML-EM iterative algorithm based on statistical model is applied to estimate Eq. (3). Taking into account the Poisson statistics in the vector X , the LS-MLEM algorithm can be written as <math display="block" id="M6"><mrow><msubsup><mi>λ</mi><mi>j</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mfrac><mrow><msubsup><mi>λ</mi><mi>j</mi><mi>k</mi></msubsup></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder></mstyle><msub><mi>h</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder></mstyle><mfrac><mrow><msub><mi>h</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>x</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>l</mi></munder></mstyle><msub><mi>h</mi><mrow><mi>i</mi><mi>l</mi></mrow></msub><msubsup><mi>λ</mi><mi>l</mi><mi>k</mi></msubsup></mrow></mfrac><mn>.</mn></mrow></math> As compared with Eq. (6), the conventional ML-EM algorithm of Eq. (1) is written as <math display="block" id="M7"><mrow><msubsup><mi>λ</mi><mi>j</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mfrac><mrow><msubsup><mi>λ</mi><mi>j</mi><mi>k</mi></msubsup></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder></mstyle><msub><mi>c</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder></mstyle><mfrac><mrow><msub><mi>c</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>p</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>l</mi></munder></mstyle><msub><mi>c</mi><mrow><mi>i</mi><mi>l</mi></mrow></msub><msubsup><mi>λ</mi><mi>l</mi><mi>k</mi></msubsup></mrow></mfrac><mn>,</mn></mrow></math> where k is the k th iteration. When the projection data includes higher-noise, the iteration for Eq. (6) is an iterating process with lower-noise, unlike Eq. (7) with higher-noise. It is well known that the ML-EM algorithm can get accurate image at lower-noise. So the noise propagated in reconstructed image, obtained using Eq. (6), can be controlled under lower-noise level. III. COMPUTER SIMULATION 1 To evaluate the properties of the LS-MLEM algorithms, we use a simulation phantom with hot and cold lesion. As shown in Fig. 1 Fig. 1 , the voxel is 3.00 mm and true image is of 121 × 121 pixels, obtained using calculation of 2420 × 2420 small pixels. The hot-to-background ratio of 5:1is assumed, and activity of the cold is zero. The radius of background is 18.15 cm, and the radius of hot and cold is 4.0 cm. The distance between hot and cold centers is 18.0 cm. A low-energy (140 keV) general-purpose parallel-beam collimator (square holes,of 2.24 mm size, 4.90 cm length, and 0.76 mm septa) is used. The SPECT head rotates around the phantom, with radius of rotation of 18.15 cm. The effects of scatter, penetration and detector response function are ignored. Using the phantom, the simulated projection data are binned into 121-element arrays for each of the 120 views over 360. Poisson noise fluctuations are then added to the simulated projection data. The projection data with SNR p = 5;10;20;30;40;50;60;70;80;90;100, i.e. the count/bin ratios of 25;100;400;900;1600;2500;3600;4900;6400;8100;10000 respectively, are used in evaluating the improved noise SNR x . In Fig. 2 Fig. 2 , the SNR x data are compared with SNR p data. The same projection data with noise are used to calculate the SNR x and SNR p . Fig. 2 Fig. 2 indicates that the SNR x is proportional to the SNR p , but the SNR x is almost 25 times larger than the SNR p . The noise of vector X is remarkably improved, and the LS-MLEM is an iteration algorithm with lower-noise. To emphasize LS-MLEM applied in projection data with higher-noise,projection data sets of SNR p = 5;10;20;30, are used in the reconstructions. The system matrix is obtained by calculating 2420 × 2420 small pixels. Both ML-EM and LS-MLEM iterative algorithms are used in reconstructing the simulated projection data. Images of 121 × 121 pixels are reconstructed from 120 projections sampled over 360. Results from the higher-noise data allow us to study properties of the iterative reconstruction methods with LS-MLEM. An important property is its convergence rate. Also, the four higher-noise data sets allow comparison of LS-MLEM with ML-EM in noise improvement. And we have a better assessment of the improved reconstructed image quality. IV. RESULT AND DISCUSSION 1 Figure 3 Figure 3 shows the reconstructed images, obtained using the ML-EM and LS-MLEM with Poisson noise added to the projection data, after the 20 th , 40 th , 90 th , 140 th , 500 th and 1000 th iteration. Figure 4 Figure 4 shows the signal-to-noise ratio as function of iteration number for the ML-EM and the LS-MLEM algorithms, plotted by the solid line and dash dot line, respectively. The signal-to-noise ratio for the full, background and hot reconstruction field are respectively plotted from first to third column. As Fig. 4 Fig. 4 shows, for ML-EM, the maximum signal-to-noise ratio in full field is at 20 th , 40 th , 90 th and 140 th iteration at SNR p = 5;10;20;30 respectively; while for LS-MLEM, the signal-to-noise increases all the way to the end. Figure 3 Figure 3 shows that the noise artifacts in the reconstructed image, obtained using ML-EM, tend to increase with iteration number. This is more obvious as the SNR p decrease which is consistent with Fig. 4 Fig. 4 . However, Figs. 3 Figs. 3 and 4 4 show that noise in the LS-MLEM-reconstructed image tends to decrease with increasing iteration number where the artifacts do not change obviously. However, it is obvious that the convergence rate of the LS-MLEM is slower than the ML-EM. To increase the convergence rate, ordered-subset is applied to LS-MLEM algorithm (called LS-OSEM in this paper). The results, obtained using LS-OSEM with 10 subsets, are plotted in Figs. 3 Figs. 3 and 4 4 . The convergence rate using LS-OSEM is improved obviously. Compared with ML-EM, the noise and artifacts from LS-OSEM do not change by much, though their increase can be seen in the reconstructed images with SNR p =5 and 10, at higher iteration number. A maximum signal-to-noise ratio is at 170 th iteration for SNR p =5. The result of SNR p = 20 and SNR p =30 consist with LS-MLEM. As shown in Figs. 3 Figs. 3 and 4 4 , the SNR and quality of image, obtained using LS-OSEM, are better than the ML-EM, for hot field especially. To further evaluate the iterative algorithms in terms of image quality, Fig. 5 Fig. 5 shows the line profile along axis of the central transverse cross section. The dotted line is of the phantom, while the other lines are defined the same as in Fig. 4 Fig. 4 . The line profile of ML-EM is at maximum signal-to-noise ratio, while the line profile of LS-MLEM and LS-OSEM are at the 1000 th and 170 th iterations, respectively. The horizontal axis is pixelnumber and vertical axis is activity. Note that Fig. 5 Fig. 5 shows the image reconstructed by ML-EM algorithm higher noise than those by the LS-MLEM and LS-OSEM algorithm, and that the line profile of LS-MLEM and LS-OSEM algorithm is the "best" estimate of assuming that the line profile satisfy the statistical model. Figs. 3 Figs. 3 , 4 4 and 5 5 also indicate that the difference between ML-EM and LS-MLEM or LS-OSEM decreases with increasing SNR p . V. CONCLUSION 1 Using standard least squares, Eq. (1), of higher-noise, can be converted to Eq. (3), of lower-noise. The SNR x is ten times more than the SNR p . This can provide an algorithm of improved noise property for the iterative image reconstruction. By incorporating ML-EM or OSEM, the LS-MLEM or LS-OSEM algorithm imposes a positivity constraint on the reconstructed image pixel values. Using the conventional OSEM algorithm for the LS-MLEM, the convergence rate is obviously improved. The LS-MLEM or LS-OSEM algorithms for SPECT imaging are advantageous in their ability of controlling higher-noise accurately to achieve improved image quality. This will not be the case when the projection data are of lower-noise, as the difference between ML-EM and LS-MLEM or LS-OSEM decreases gradually with the noise, however, the projection noise is heavily contaminated by statistical noise, practically. The LS-MLEM or LS-OSEM is meaningful in SPECT application. Further efforts will be made in modeling the imaging process for compensating image degrading factors, such as attenuation, scatter and spatial resolution. References Shepp L A and Vardi Y . IEEE T Med Imaging , 1982 , 2 : 113 - 122 . IEEE T Med Imaging MHudson H and Larkin R S . IEEE T Med Imaging , 1994 , 4 : 601 - 609 . IEEE T Med Imaging Huesman R H . Phys Med Biol , 1984 , 5 : 543 - 552 . Phys Med Biol Llacer J , Veklemv E , Coakley K I , et al . IEEE T Med Imaging , 1993 , 2 : 215 - 231 . IEEE T Med Imaging Green P J . IEEE T Med Imaging , 1990 , 1 : 84 - 93 . IEEE T Med Imaging Tsui B M W , Zhao X D , Frey E C , et al . IEEE T Nucl Sci , 1991 , 6 : 1766 - 1772 . IEEE T Nucl Sci La V and Grangeat P . Phys Med Biol , 1998 , 43 : 715 - 727 . Phys Med Biol Zeng G L , Gagnon D , Natterer F , et al . IEEE T Nucl Sci , 2003 , 5 : 1590 - 1594 . IEEE T Nucl Sci 10.13538/j.1001-8042/nst.25.030302 Supported by the Priority Academic Program Development of Jiangsu College Education 2013-10-10 2013-12-18 2014-06-20 2014-06-20 2021-04-11 2014 Nuclear Science and Techniques 3 25