I. INTRODUCTION

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, 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].

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, 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-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 solutions 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

In SPECT, the process of obtaining projection data is as follows:

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 cil representing the relative contribution of l^th voxel to projection bin pi. 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

The squared residual error E(λ) shall be minimized. According to standard least squares, we have

where C^T is transpose matrix of the system matrix C. Defining X = C^TP, its elements xj are determined by

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 xj in Eq. (4) is a linear superposition of the projection data with Poisson noise. Therefore, the vector X, with its elements xj, 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/2SNR_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^TC be a constant matrix based on system matrix, its elements hjl is determined by

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

As compared with Eq. (6), the conventional ML-EM algorithm of Eq. (1) is written as

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

To evaluate the properties of the LS-MLEM algorithms, we use a simulation phantom with hot and cold lesion. As shown in 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, 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 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.

Fig. 1.Full-screen View

The true phantom.

Fig. 2.Full-screen View

SNR_x versus SNR_p.

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

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.

Fig. 3.Full-screen View

Reconstructed images of the phantom in Fig. 1 (The LS-OSEM results will be discussed later).

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 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.

Fig. 4.Full-screen View

SNR versus iteration number (—, ML-EM; -⋅-, LS-MLEM; ---, LS-OSEM).

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. However, Figs. 3 and 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 and 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 and 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 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. 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 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, 4 and 5 also indicate that the difference between ML-EM and LS-MLEM or LS-OSEM decreases with increasing SNR_p.

Fig. 5.Full-screen View

Profile (⋯, true phantom, —, ML-EM; -⋅-, LS-MLEM; ---, LS-OSEM).

V. CONCLUSION

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.