Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Med Biol. Author manuscript; available in PMC 2013 June 24.
Published in final edited form as:
PMCID: PMC3690946

Quantitative evaluation study of four-dimensional gated cardiac SPECT reconstruction


In practice gated cardiac SPECT images suffer from a number of degrading factors, including distance-dependent blur, attenuation, scatter, and increased noise due to gating. Recently we proposed a motion-compensated approach for four-dimensional (4D) reconstruction for gated cardiac SPECT, and demonstrated that use of motion-compensated temporal smoothing could be effective for suppressing the increased noise due to lowered counts in individual gates. In this work we further develop this motion-compensated 4D approach by also taking into account attenuation and scatter in the reconstruction process, which are two major degrading factors in SPECT data. In our experiments we conducted a thorough quantitative evaluation of the proposed 4D method using Monte Carlo simulated SPECT imaging based on the 4D NURBS-based cardiac-torso (NCAT) phantom. In particular we evaluated the accuracy of the reconstructed left ventricular myocardium using a number of quantitative measures including regional bias-variance analyses and wall intensity uniformity. The quantitative results demonstrate that use of motion-compensated 4D reconstruction can improve the accuracy of the reconstructed myocardium, which in turn can improve the detectability of perfusion defects. Moreover, our results reveal that while traditional spatial smoothing could be beneficial, its merit would become diminished with the use of motion-compensated temporal regularization. As a preliminary demonstration, we also tested our 4D approach on patient data. The reconstructed images from both simulated and patient data demonstrated that our 4D method can improve the definition of the LV wall.

Keywords: 4D reconstruction, gated cardiac SPECT, motion-compensated temporal smoothing, attenuation and scatter compensation

1. Introduction

Single photon emission computed tomography (SPECT) is an important diagnostic imaging technique in use for diagnosis and evaluation of cardiac diseases. In gated cardiac SPECT, the data acquisition is synchronized to the electrocardiogram (ECG) signal, and a sequence of 3D images is reconstructed instead of a single static image. Consequently, gated SPECT can provide valuable ventricular function information such as ejection fraction, wall motion, and wall thickening (Garcia 1996). However, the effectiveness of SPECT often suffers from a number of degrading factors, ranging from decreased image contrast caused by scatter, image artifacts due to attenuation, poor spatial resolution due to distance-dependent blur, to increased noise in gated images.

In recent years there has been considerable interest in development of image processing methods for noise reduction in gated cardiac imaging (King & Miller 1985, Narayanan et al. 2000, Peter et al. 2001, Brankov et al. 2005, Klein et al. 1997, Lalush & Tsui 1998, Gilland et al. 2002, Cao et al. 2003, Gravier et al. 2006). Generally speaking, these methods can be divided into two broad categories: 1) pre- or post-filtering methods, and 2) joint reconstruction methods. Both approaches aim to exploit the statistical correlation (i.e., similarity) among the different gated frames in a cardiac sequence. In the pre- or post-filtering approach, the individual gated frames are reconstructed independently, either preceded by or followed by filtering along the temporal direction (e.g., (King & Miller 1985, Narayanan et al. 2000, Peter et al. 2001, Brankov et al. 2005)). In the joint reconstruction approach, the gated frames are reconstructed jointly from the entire sequence of data (e.g., (Klein et al. 1997, Lalush & Tsui 1998, Gilland et al. 2002, Cao et al. 2003, Gravier et al. 2006)).

In our recent work (Gravier et al. 2006), we developed a spatio-temporal reconstruction method in which estimated cardiac motion was used to enforce the correlation among the gated frames. We demonstrated the feasibility of this approach based on several quantitative measures using the 4D gated mathematical cardiac-torso (gMCAT) phantom (Pretorius et al. 1999). Our results showed that incorporation of motion-compensated temporal regularization yields effective noise reduction, resulting in improved perfusion-defect detection in the reconstructed images while avoiding significant cross-frame blurring.

Encouraged by this success, we now further develop and evaluate this approach in a more realistic setting, which is a necessary step toward our ultimate goal of producing a useful clinical tool. In our previous work (Gravier et al. 2006), we excluded attenuation and scatter so we could isolate and understand the effect of motion compensation by itself. Of course, to be realistic, the algorithm must account for attenuation and scatter. In this paper we now include both, and perform extensive evaluations of the algorithm in this more-realistic situation, including both simulated and patient data.

In SPECT imaging, both attenuation and scatter can have severe effects on the quality of the image data (Jaszczak et al. 1981, Zaidi & Hasegawa 2003, deVries & King 1994, Hademenos et al. 1993). If not properly corrected, both factors are known to cause artifacts in the reconstruction, which could exhibit as decreased contrast (Heller et al. 1997) or even spurious perfusion defects (Germano et al. 1994).

In this study, we account for the attenuation factor in the reconstruction algorithm by modeling the attenuation map in the system matrix. For scatter compensation, the scatter component is modeled in the likelihood function of the acquired data (as indicated later in Eq. (1) and thereafter). The resulting 4D reconstruction algorithm will now take into account the following four major degrading factors in gated SPECT: low sensitivity (more severe due to gating), distance-dependent spatial resolution, attenuation, and scatter. This allows us to evaluate the benefit of our 4D approach in a more realistic setting.

Moreover, in our evaluation studies, Monte Carlo simulation (Ljungberg & Srand 1989) is used to simulate the acquisition data using the more recently developed 4D NURBS-based cardiac-torso (NCAT) phantom (Segars 2001). Owing to both the complex photon processes inside the phantom and the interactions between the photons and the gamma camera simulated by the Monte Carlo method, this allows for an unbiased evaluation of the reconstruction performance in the presence of a potential mismatch between the actual imaging process and the analytic imaging model used in the reconstruction algorithm.

For the numerical reconstruction we now apply a modified block sequential regularized expectation-maximization (BSREM) algorithm (Ahn & Fessler 2003), which is globally convergent (Ahn & Fessler 2003) and much faster than the one-step late (OSL) algorithm used in (Gravier et al. 2006).

In our evaluation studies, we use the following three quantitative criteria: 1) reconstruction accuracy of the myocardium, quantified by mean squared error (MSE) and bias-variance plot; 2) cardiac perfusion defect detection, measured by a channelized Hotelling observer (CHO) (Myers & Barett 1987); and 3) time activity curves (TACs) and ejection fraction (EF) of the left ventricle, a clinical measure of left ventricular function. As a preliminary demonstration of the proposed approach, we also provide reconstruction results from a clinical acquisition.

The rest of the paper is organized as follows. The imaging model for gated SPECT and our proposed reconstruction method are described in Sect. II. The evaluation methods are given in Sect. III, and the evaluation results are presented in Sect. IV. Conclusions are given in Sect. V.

2. Reconstruction of gated cardiac SPECT

2.1. Imaging model

In gated cardiac SPECT, the acquired projection data are binned into K cardiac-gate intervals by using the ECG signal. The imaging data are described by the following model


in which gk, fk and rk are vectors representing the acquired data, original image, and expected scatter component, respectively, of the kth gate interval, E[·] is the expectation operator, and H is a matrix describing the imaging system called the system matrix. In H each element hij represents the probability that a photon emitted at voxel location j is detected at detector bin i without being scattered.

The goal is to estimate the images fk from the data gk, k = 1, 2,…, K. Owing to the ill-conditioned nature of the system matrix H, a direct inversion of (1) is highly sensitive to the imaging noise. In 4D reconstruction, we seek a simultaneous solution of all the frames fk from (1), aiming to compensate for H by taking into account both the spatial and temporal correlations among the gate frames.

In this study the elements hij of the system matrix H in (1) are modeled including both the distance-dependent point spread function (PSF) and the attenuation effect. Specifically, let hij0 denote the elements of the system matrix with distance dependent spatial resolution, but without attenuation. Let μ(x, y, z) denote the attenuation map of the object. Then the attenuation effect is modeled in the system matrix as


where the line integral is carried out along the line segment between image voxel j and detector bin i. In our experiments, the attenuation map from NCAT phantom was used for the simulation data; for the patient data μ(x, y, z) was estimated from the transmission images.

It is noted that in our experiments the projection data were generated by Monte Carlo simulation using the SIMIND program (Ljungberg & Srand 1989) rather than by an analytic projection model. This allows for a more accurate modeling of the SPECT imaging process (as reflected by (1)) in our evaluation studies. For reconstruction, the analytic system matrix H was used. This provides a more objective and realistic evaluation of the reconstruction algorithms.

2.2. Maximum a posteriori (MAP) estimate

For convenience, let F[equivalent](f1T,f2T,,fKT)T, i.e. a vector consisting of all the gate frames; similarly, let G[equivalent](g1T,g2T,,gKT)T. Then, the maximum a posteriori (MAP) estimate of F is obtained as:


where p(G; F) is the likelihood function of G parameterized by F, and p(F) is a prior distribution on F. Below we explain the terms p(G; F) and p(F) separately.

We assume a Poisson likelihood function for the noise model in (1) for SPECT. Then we have


where (Hfk)(i), gk(i) and rk(i) denote the ith elements of (Hfk), gk and rk, respectively, and M is the total number of detector bins.

As mentioned above, the matrix H in (4) is modeled to include both the distance-dependent PSF and the attenuation effect of a SPECT system. Moreover, the scatter component rk has been directly taken into account in the likelihood function as in (Lange & Carson 1984). Such a formulation will lead to automatic compensation for these degrading factors in the resulting iterative reconstruction algorithm.

The prior distribution p(F) is used to ameliorate the noise effect in the reconstructed images. Traditionally, p(F) is defined to enforce local spatial smoothness. In gated SPECT, the strong temporal correlation among the gated frames also warrants temporal regularization. Here we follow closely our earlier approach in (Gravier et al. 2006). Specifically, we assume a Gibbs prior for p(F) of the following form

p(F)[proportional, variant]exp[βsUs(F)βtUt(F)]

where Us(F), Ut(F) are two energy functions defined in space and time, respectively, and βs, βt are their corresponding scalar weighting parameters. The spatial term Us(F) is used to exploit the similarity among neighboring voxels (Gravier et al. 2006). The gate term Um(F) is used to exploit the similarity among different gate intervals. Specifically, it is defined to enforce the consistency across the different gate frames by following the motion trajectories of the image voxels (Gravier et al. 2006, Jin et al. 2007). The same definitions as in (Jin et al. 2007) are used for these prior terms in this study, which are omitted here for brevity.

As to be illustrated later in Sect. 3.3, by properly setting of the parameters βs and βt in the prior term p(F), we can obtain reconstruction algorithms corresponding to several different types of spatio-temporal regularization.

With all the terms defined, the MAP estimate in (3) can then be solved by minimizing the following objective function


2.3. Reconstruction algorithm

For the minimization problem in (6), the modified block sequential regularized expectation-maximization (BSREM) algorithm (Ahn & Fessler 2003) is adapted in this study. In our previous work (Gravier et al. 2006), the one-step late (OSL) algorithm was used. The BSREM algorithm is known to be globally-convergent and typically faster than non-ordered-subset algorithms. As in a standard ordered-subset algorithm (Herman & Meyer 1993), the set of projection data G is divided into a number of subsets Sd, d = 1,…,D, which results in J(F)=d=1DJd(F). Then the image frames F are updated as follows


where n is the iteration index, d is the sub-iteration index, p is a projection operator onto θ [equivalent] {F : [sm epsilon]FU[sm epsilon]} with [sm epsilon] being a small positive number and U being an upper bound defined as in (Ahn & Fessler 2003), νn is a relaxation parameter, [down-pointing small open triangle] is the gradient operator, and finally, Λ(Fn,d−1) is a diagonal scaling matrix defined below.

In (7), the diagonal entry of Λ(Fn,d−1) for the update on fk(j) is given by:


where pj is chosen as


In our evaluation study, D = 16 was used; the subsets were formed such that two consecutively updated subsets were separated as far as possible. The number of iterations of BSREM used in this study was 10.

3. Evaluation methods

3.1. Data sets

3.1.1. Phantom simulation data

In our simulation study the 4D NURBS-based cardiac-torso (NCAT) 2.0 phantom (Segars 2001) was used to generate the source and attenuation distribution. A transmural perfusion defect with circumferential width 45° and long-axis width about 25mm was introduced in the septal wall of the left ventricle (LV). The intensity reduction of the defect is 20%. The SIMIND Monte Carlo package (Ljungberg & Srand 1989) was used to simulate gated SPECT imaging with Tc99m labeled sestamibi as the imaging agent. The SPECT system simulated was the Picker Prism3000 with a low-energy high-resolution (LEHR) collimator. The projections were 64 by 64 bins with a pixel size of 0.634 cm. For a circular camera rotation of 28.5 cm radius, 64 projection sets were collected for each gate frame. A total of 16 gates were used. Two energy windows were used in SIMIND as in (Ogawa et al. 1994): the photopeak window was 20% (28keV) centered at 140keV, and a 3.5 keV window abutted to the lower side of the photopeak window. Poisson noise was introduced at a level of 8 million total acquired counts as in a typical clinical acquisition. The scattered counts were about 34% of the total counts.

In SIMIND the Monte Carlo method was used to simulate the complex photon processes in both the phantom and the interactions of photons with the gamma camera. This would remove any potential bias that may arise when the same analytical projection model is used for both data generation and reconstruction.

Finally, as the ground truth for comparison, the gated volumetric images were reconstructed using the ML-OSEM algorithm (10 iterations with 16 subsets) from simulated noiseless projection data of the NCAT without attenuation and scatter (by turning off attenuation and scatter in SIMIND); depth-dependent blurring was included both in the simulated data and in the reconstruction algorithm. To minimize the noise effect, we used a total of 640 million events in the Monte Carlo simulation, from which the ideal reconstruction was obtained. This represents the ideal case of perfect acquisition where none of attenuation, scattering, and Poisson noise are involved in the sinogram data (denoted as “Ideal” below for reference). Thus, it represents what would be achievable in the scenario of both perfect attention and scatter correction and perfect noise filtering. By comparing to this noiseless reconstruction, we will learn how far the proposed method is from this ideal situation. These images were used as reference for quantitative evaluation of reconstructed images as described below.

3.1.2. Clinical data

As a preliminary demonstration, we also used a set of clinical data. This data set was acquired from a female patient by an IRIX system (Gagnon et al. 1999) with 68 projections (3-degree steps) and a 128 × 128 matrix. The pixel size was 0.467 cm. The acquisition started from right anterior oblique, passed through anterior and left anterior oblique, and ended at left lateral oblique. A total of 8 gates were used. The photopeak window was 15% centered at 147.5 keV, and a second 5% window centered at 133.2 keV was used. The total number of counts acquired was 23.8 millions. For attenuation correction, the attenuation map was estimated from the transmission images of the patient. Note that in this set of clinical data the photopeak window was set to be centered at 147.5 keV so that it was asymmetric to the high side in order to decrease scatter within the acquired photopeak window data (Atkins et al. 1968).

3.2. Evaluation criteria

We evaluated the different reconstruction methods using several quantitative criteria. They are organized into the following three categories: 1) reconstruction accuracy measures of the left ventricle (LV) myocardium, 2) detectability of cardiac perfusion defects using a channelized Hotelling observer (CHO) (Myers & Barett 1987); and 3) LV function measures. We describe these measures next.

3.2.1. Reconstruction accuracy of the LV myocardium

To quantify the overall accuracy of the reconstructed LV myocardium, we computed the mean squared error (MSE) of a 30 × 28 × 20 volumetric region containing the entire LV, of which a 2D slice is shown in Fig. 1. The MSE was computed for all 16 gates. In addition, to quantify the reconstruction accuracy of the LV wall, we conducted bias-variance analyses of two selected regions of interest (ROI) (one normal and one defect, as indicated in Fig. 1). The ideal images described above were used as reference in these quantitative measures.

Figure 1
Regions of interest (ROIs) used for different quantitative measurements: normal and defect ROIs in bias-standard deviation plot, and time activity curve. The center location for the CHO detector is indicated by “*”.

We also quantified the uniformity of the LV wall based on a short-axis view. Specifically, a set of 32 voxels were selected along the heart wall in one short-axis slice; the mean of the intensity values at this set of voxels was computed, and their average deviation from the mean was then computed as a measure of the wall uniformity. Specifically, the wall uniformity is defined as

U=i=1M[mid ]fROI(i)fROI[mid ]MfROI×100%

where {fROI(i), i = 1,…, M} denotes the selected set of voxels on the heart wall, and fROI is the mean value of this set of voxels.

3.2.2. Cardiac perfusion defect detection

We applied a CHO to quantify the detectability of the introduced perfusion defect in the reconstructed images. In our implementation, a 28 × 28 region centered around the lesion region (as shown in Fig. 1) was bilinearly interpolated into a 140 × 140 image. We used four rotationally symmetric, non-overlapping input channels for the CHO with the following passbands: (0.03125, 0.0625), (0.0625, 0.125), (0.125, 0.25) and (0.25, 0.5) cycles/pixel. These channels were similar to those used in (Myers & Barett 1987). Internal noise was introduced as in (Narayanan et al. 2002). The detection performance was summarized using the area under the receiver operating characteristic (ROC) curve (Metz 1986). For each method, a total of 60 noise realizations of reconstructed images were used: 30 with the defect present and 30 with the defect absent; the statistics were then estimated for the CHO for obtaining the area under the ROC curve as in (Narayanan et al. 2002). These ROC studies were “signal-known exactly” and “background-known exactly” observer studies (Narayanan et al. 2002). In our experiments, the CHO was computed on the first gate (near end-diastole) interval.

3.2.3. LV function

We also quantified the reconstruction using ejection fraction (EF), a measure pertinent to clinical assessment of LV functions. This parameter was computed using the clinical software package 4DM-SPECT.

It is noted that with the clinical software the image analysis parameters may not be optimized for the proposed reconstruction methods, owing to the different resolution properties of the images. To accommodate this, we applied a training step (with a different set of noise realizations) prior to using the software, based on which the reconstructed images were first pre-processed with a lowpass Gaussian filter (std=1).

3.2.4. Myocardial TAC results

In addition, to demonstrate the effect of potential temporal smoothing on cardiac motion, we also calculated the time activity curves (TACs) of an ROI at the base of the LV myocardium shown in Fig. 1. As the wall moves in and out of this ROI during the beating cycle, its average intensity will vary accordingly, and thus, serves as a good indicator on the degree on temporal smoothing.

The normalized cross-correlation coefficients (CC) between a reconstructed TAC and its ideal reconstructed one are also computed as follows

CC(y,y^)=(yy)T(y^y^¯)yy[center dot]y^y^¯

where y = [y(1),…, y(K)]T and y^=[y^(1),,y^(K)]T are the ideal and reconstructed TACs on the selected ROI; y and y^¯ are the mean values of y and y^.

3.3. Reconstruction methods

We tested the reconstruction algorithms using different settings of βs and βt. Specifically, 1) βs ≠ 0 and βt = 0, i.e., MAP reconstruction with only spatial smoothing (MAP-S); 2) βs = 0 and βt ≠ 0, i.e., MAP reconstruction with only temporal smoothing; and 3) βs ≠ 0 and βt ≠ 0, i.e., MAP reconstruction with both spatial and temporal smoothing. Note that when βs = 0 and βt = 0, it amounts to reconstruction with both scatter and attenuation corrections, but no prior. This corresponds to maximum likelihood (ML) reconstruction.

The filtered TEW method (Ogawa et al. 1991, King et al. 1997) was used to estimate the scatter component in the projection data (i.e. rk in (6)). For the NCAT data, an estimate of the scatter counts within the photopeak window was first obtained by multiplying the counts in the 3.5 keV window by 4. This estimate was further filtered with an order-3 2D Butterworth lowpass filter with a cutoff frequency of 0.2 cycles/cm as suggested in (King et al. 1997). For attenuation correction, we used the attenuation map generated by the NCAT program. For the clinical data, the scatter component was estimated similarly except that the estimate of the scatter counts within the photopeak window was obtained by multiplying the counts in the secondary window by 1.66.

To determine the motion prior, we first applied the BSREM reconstruction algorithm with only spatial smoothing (MAP-S) to get an initial estimate of the gate frames. Afterwards a 3D version of the optical flow method (Horn & Schunck 1981) was applied to estimate the motion between the different gate frames as in our previous work (Gravier et al. 2006).

For comparison, we also tested a clinical spatio-temporal processing method. In this method, the images were first reconstructed by using filtered backprojection (FBP), then processed spatially with a Butterworth filter of order 2.4 and cutoff frequency of 0.2 cycles/voxel, and further processed with a temporal filter with impulse response 1/4, 2/4, 1/4. We will refer to it as the ST121 method below. To compensate for the lack of attenuation and scatter correction, in our experiments the ST121 reconstruction was normalized such that it has the same total intensity in the reconstructed volume as the ideal images.

4. Results and discussions

4.1. Phantom simulation data

4.1.1. Reconstruction accuracy of the myocardium

The MSE results of the reconstructed LV myocardium are summarized in Fig. 2 for the different methods. These results were averaged over 30 noise realizations. In Fig. 2, each curve was obtained by varying the spatial parameter βs ([0, 0.5, 2, 5, 10]×10−4) while the temporal parameter βt was held at a constant value. For comparison, ST121 is also shown.

Figure 2
Mean square error (MSE) of the reconstructed LV myocardium with different settings of βs and βt. These results were obtained from 30 noise realizations. For comparison the MSE obtained by ST121 reconstruction is also shown.

The results in Fig. 2 show that motion-compensated 4D reconstruction (βt ≠ 0) was more accurate than spatial only MAP reconstruction (βt = 0), and both were more accurate than ST121. The best MSE value (10.67) was obtained with βs = 0 and βt = 7 × 10−4, nearly 8 times smaller than that of ST121.

Moreover, when there was little or no temporal smoothing (e.g., βt = 1 × 10−4), additional spatial smoothing could resulted in further improvement in reconstruction accuracy. However, the merit of spatial smoothing seems to diminish with increased temporal smoothing (e.g., βt = 7×10−4). This can be explained as follows. The LV wall is only a few voxels thick in some gates; thus, while spatial smoothing can be effective for noise reduction, it can also lead to spatial blurring of the LV wall (as evidenced from reconstructed images shown later in Figs. 6 and and7).7). With motion-compensated temporal smoothing, the benefit of spatial smoothing is outweighed by its deleterious impact on resolution of the heart wall.

Figure 6
Reconstructed images of a transverse slice of four gates (#1, #5, #9 and #13) by different reconstruction methods: ST121, ML, MAP-S (βs = 5 × 10−4), 4D reconstruction with βs = 5 × 10−5, βt = 7 × ...
Figure 7
Reconstructed images of four gates (#1, #5, #9 and #13) by different reconstruction methods in short-axis view: ST121, ML, MAP-S (βs = 5 × 10−4), 4D reconstruction with βs = 5 × 10−5, βt = 7 × ...

Next, we show in Fig. 3 plots of bias vs. standard deviation of the two ROIs in Fig. 1. A total of 30 noise realizations were used; the mean intensity of each ROI was computed for each realization, and its bias and variance values were obtained from the multiple noise realizations. As earlier in Fig. 2, each curve in Fig. 3 was obtained by varying the spatial smoothness parameter βs ([0, 0.5, 2, 5, 10]×10−4) while the temporal parameter βt was held at a constant value.

Figure 3
Bias-standard deviation plots for the normal and defect ROIs reconstructed with different settings of βs and βt. These results were obtained from 30 noise realizations. For reference the result is also shown for ST121.

In Fig. 3, each curve shows a decreasing trend in the ordinate (standard deviation) and an increasing trend in the abscissa (bias) as βs is increased. Moreover, use of increased temporal smoothing (βt > 0) could further decrease both the variance and bias, except for only a slight increase in bias when βs = 0 and βt increased from 3×10−4 to 7 × 10−4 for the normal ROI. The ST121 method could achieve a small variance but with a rather large bias.

We note that the MSE was used to measure the error in the whole myocardial volume, whereas the bias-variance was used to quantify the accuracy of the selected ROIs in the heart wall. Therefore, it is not surprising that these measures would lead to slightly different optimal values. Nevertheless, the above results demonstrate that the performance reaches a plateau for βt in the range from 3 × 10−4 to 2 × 10−3.

Finally, we show in Table 1 the uniformity results of the LV wall obtained by the different methods. These results were computed from the reconstructed images of 30 noise realizations by each method (the standard deviation values are given in parentheses). These results show that motion-compensated 4D reconstruction could also lead to more accurate preservation of the uniformity of the LV wall in the reconstructed images (as evidenced from reconstructed images shown later in Figs. 6, ,77 and and88).

Figure 8
Reconstructed images of four gates (#1, #5, #9 and #13) by different reconstruction methods in short-axis view (with defect): ST121, ML, MAP-S (βs = 5 × 10−4), 4D reconstruction with βs = 5 × 10−5, β ...
Table 1
Uniformity measure of the LV wall reconstructed different methods: 1) spatio-temporal smoothing (ST121), 2) maximum likelihood (ML), 3) MAP with spatial only smoothing (MAP-S), and 4) 4D with βs = 5 × 10−5, βt = 3 × ...

4.1.2. Cardiac perfusion defect detection

We summarize the lesion detectability CHO results for the different reconstruction algorithms in Fig. 4, where the area under the ROC curve (Az) is plotted for different parametric settings. Clearly, motion-compensated 4D reconstruction (βt ≠ 0) achieved the best detection performance, outperforming both spatial only MAP reconstruction (βt = 0) and ST121. Interestingly, these results are also consistent with the MSE and bias-standard deviation results presented above; in particular, use of spatial smoothing could improve the lesion detectability when there was no or little temporal smoothing (e.g., βt = 1 × 10−4), however, its merit would diminish with increased temporal smoothing (e.g., βt = 7 × 10−4).

Figure 4
CHO detection Az results of reconstruction with different settings of βs and βt.

The above results demonstrate that 4D reconstruction could potentially lead to improved detectability of perfusion defects in individual gated frames. In current clinical practice, however, perfusion defects are evaluated primarily using ungated images because of the high noise level in conventional gated images. With noise suppression, 4D reconstruction may open the door to study the feasibility of using gated images for defect detection, which may enable higher sensitivity since it would avoid blurring of the defect due to cardiac motion. This is certainly an interesting topic worthy of further clinical investigation in our future research.

4.1.3. LV function

The EF results calculated from images reconstructed by different methods are summarized in Table 2. These results were averaged from 30 noise realizations (with standard deviations given in parentheses). As reference, the ground truth of these parameters of the NCAT phantom is also listed. It is observed from these results that all the different methods could lead to fairly accurate estimates of the EF (with the mean all within 5% of the true EF value); the 4D results also show a smaller variance, thanks to reduced noise in the reconstructed images.

Table 2
LV ejection fraction (EF) values from different methods: 1) spatio-temporal smoothing (ST121), 2) MAP with spatial only smoothing (MAP-S), and 3) 4D with βs = 5 × 10−5, βt = 3 × 10−4. These results were ...

4.1.4. Myocardial TAC results

In Fig. 5 we show the TACs of the ROI shown in Fig. 1 across the different gates from 30 noise realizations, where the average TACs are shown along with their standard deviations for the different reconstruction methods; for comparison purposes, in this plot the ST121 images were normalized such that the myocardium intensity (average from 30 noise realizations) has the same maximum as that of the ideal reference. For clarity, 4D results are shown only for the case of βs = 5 × 10−5 and βt = 3 × 10−4. As can be seen, the TACs of 4D exhibit less bias when compared to that of ST121, particularly at the end-systole (ES) stage (6-8 gates). This is consistent with the bias-standard deviation results shown earlier in Fig. 3. We also computed the normalized cross-correlation coefficients (CC) between the reconstructed and the ideal TACs, which are given in Table 3. These results were based on 30 noise realizations.

Figure 5
Time activity curves across 16 gates obtained with different methods: clinical spatio-temporal smoothing (ST121); maximum likelihood (ML); 4D reconstruction with βs = 5 × 10−5, βt = 3 × 10−4. A total of ...
Table 3
Normalized cross-correlation coefficients of reconstructed TACs by different methods: 1) spatio-temporal smoothing (ST121), 2) maximum likelihood (ML), 3) MAP with spatial only smoothing (MAP-S), and 4) 4D with βs = 5 × 10−5, ...

4.1.5. Reconstructed images

In Fig. 6 we show a set of typical images reconstructed by the different methods, namely, Ideal, ST121, MAP-S and 4D with βs = 5 × 10−5 and βt = 7 × 10−4; for clarity, only one transverse slice of the LV myocardium (slice #36) is shown for four selected gates. Similarly, in Fig. 7 we show a short-axis slice of the reconstructed images by the different methods. In addition, in Fig. 8 we show a set of reconstructed images when the perfusion defect was present.

From these results it can be seen that 4D images suffer from far less blur than that of ST121. The LV wall in 4D is notably better defined, particularly near the ES stage (gate #9). The images in Fig. 6 show that the septal wall of the LV is also more accurately reconstructed in 4D; as a consequence, the short-axis view images in Fig. 7 show that the image intensity is more uniform along the LV wall for 4D. This is consistent with the results shown earlier in Table 1.

Finally, we note that in this study the system matrix H in the imaging model in (1) is treated as constant, which is consistent with current clinical practice where a single attenuation map is acquired using transmission tomography. A more accurate model would be to also include potential variations associated with each cardiac phase. This could be of interest for further investigation.

4.2. Clinical data reconstruction

In Fig. 9 we show a set of reconstructed images from the patient data by our 4D method and ST121. Owing to the lack of ground truth, the parameters for 4D were chosen empirically (βs = 3 × 10−5 and βt = 3 × 10−5); for comparison, results were also shown for iterative reconstruction with no prior (i.e.,βs = 0 and βt = 0). As can be seen, the LV wall is better defined in 4D, and appears to be more uniform. In particular, the heart wall is better resolved from the high bowel activities nearby.

Figure 9
Reconstructed images of four gates (#1, #3, #5 and #7) from a set of patient data with ST121, ML and 4D (βs = 3 × 10−5, βt = 3 × 10−5).

5. Conclusions

We developed a motion-compensated 4D reconstruction approach for gated cardiac SPECT, aiming to improve the accuracy of reconstruction in the presence of degradation factors such as increased noise due to gating, distance-dependent blur, attenuation and scatter. With all these factors considered, the gate frames in a sequence were jointly reconstructed by using a modified block sequential regularized expectation-maximization (BSREM) algorithm. The proposed algorithm was quantitatively evaluated by using several criteria based on the accuracy of the LV myocardium. Our results demonstrated that motion-compensated 4D reconstruction can lead to improved accuracy of the myocardium. Interestingly, the results showed that while spatial smoothing could lead to improvement in reconstruction accuracy when there was little or no temporal smoothing, its benefit would diminish with increased motion-compensated temporal smoothing. The reconstructed images from both simulated and patient data demonstrated that our 4D method could improve the definition of the LV wall (i.e., less distortion and artifacts). Encouraged by these promising results, we plan to conduct a subsequent clinical evaluation using clinical acquisitions and expert observers.


This work was supported in part by the National Institutes of Health under grant HL65425.


  • Ahn S, Fessler JA. Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms. IEEE Transactions on Medical Imaging. 2003;22:613–626. [PubMed]
  • Atkins FB, Beck RN, Hoffer PN, Palmer D. Dependence of optimum baseline setting on scatter fraction and detector response function. Medical Radioisotope Scintigraphy, IAEA. 1968:101–118.
  • Brankov JG, Yang Y, Wernick MN. Spatio-temporal processing of gated cardiac spect images using deformable mesh modeling. Med. Phys. 2005;32(9):2839–2849. [PubMed]
  • Cao Z, Gilland DR, Mair BA, Jaszczak RJ. Three-dimensional motion estimation with image reconstruction for gated cardiac ECT. IEEE Transactions on Nuclear Science. 2003;50(3):384–388.
  • deVries DJ, King MA. Window selection for dual photopeak window scatter correction in Tc-99m imaging. IEEE Transaction on Nuclear Science. 1994;41:2771–2778.
  • Gagnon D, Tung CH, Zeng L, Hawkins WG. Design and early testing of a new medium-energy transmission device for attenuation correction in SPECT and PET. Proc. of 1998 IEEE Medical Imaging Conderence.1999. pp. 1349–1353.
  • Garcia EG. Imaging guidelines for nuclear cardiology procedures part I. J. Nucl. Card. 1996;3:G1–G46.
  • Germano G, Chua T, Kiat H, Areeda JS, Berman DS. A quantitative phantom analysis of artifacts due to hepatic activity in technetium-99m myocardial perfusion SPECT studies. J. Nucl. Med. 1994;35(2):356–359. [PubMed]
  • Gilland DR, Mair BA, Bowsher JE, Jaszczak RJ. Simultaneous reconstruction and motion estimation for gated cardiac SPECT. IEEE Transactions on Nuclear Science. 2002;49:2344–2349.
  • Gravier E, Yang Y, King MA, Jin M. Fully 4d motion-compensated reconstruction of cardiac SPECT images. Phys. Med. Biol. 2006;51:4603–4619. [PubMed]
  • Hademenos GJ, King MA, Ljungberg M, Zubal IG, Harrell CR. A scatter correction method for Tl-201 images: a Monte Carlo investigation. IEEE Transaction on Nuclear Science. 1993;40:1179–1186.
  • Heller EN, DeMan P, Liu YH, Dione DP, Zubal IG, Wackers FJ, Sinusas AJ. Extracardiac activity complicates quantitative cardiac SPECT imaging using a simultaneous transmission-emission approach. J. Nucl. Med. 1997;38:1882–1890. [PubMed]
  • Herman GT, Meyer LB. Algebraic reconstruction techniques can be made computationally efficient. IEEE Trans. Med. Imag. 1993;12:600–609. [PubMed]
  • Horn BKP, Schunck BG. Determining optical flow. Artificial Intelligence. 1981;17:185–203.
  • Jaszczak RJ, Coleman RE, Whitehead FR. Physical factors affecting quantitative measurements using camera-based single photon emission computed tomography (SPECT) IEEE Transaction on Nuclear Science. 1981;28:69–80.
  • Jin M, Yang Y, Brankov JG, Wernick MN, Feng B, King MA. Four-dimensional reconstruction of gated cardiac SPECT with attenuation and scatter compensation. Proc. IEEE Intr. Symp. Biomed. Imag.:Macro to Nano. 2007:169–172.
  • King MA, deVries DJ, Pan TS, Pretorius PH, Case JA. An investigation of the filtering of TEW scatter estimates used to compensate for scatter with ordered subset reconstruction. IEEE Transactions on Nuclear Science. 1997;44:1140–1145.
  • King MA, Miller TR. Use of nonstationary Wiener filter in nuclear medicine. Eur. J. Nucl. Med. 1985;10:458–461. [PubMed]
  • Klein GJ, Reutter BW, Huesman RH. Non-rigid summing of gated PET via optical flow. IEEE Transactions on Nuclear Science. 1997;44:1509–1512.
  • Lalush DS, Tsui BMW. Block-iterative techniques for fast 4D reconstruction using a prior motion models in gated cardiac SPECT. Phys. Med. Biol. 1998;43:875–886. [PubMed]
  • Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J. Comp. Assist. Tomogr. 1984;8:306–318. [PubMed]
  • Ljungberg M, Srand SV. A Monte Carlo program simulating scintillation camera imaging. Comp. Meth. Prog. Biomed. 1989;29:257–272. [PubMed]
  • Metz CE. ROC methodology in radiologic imaging. Investigative Radiology. 1986;21:720–733. [PubMed]
  • Myers KJ, Barett HH. Addition of a channel mechanism to the ideal-observer model. Journal of the Optical Society of America A. 1987;4(12):2447–2457. [PubMed]
  • Narayanan MV, Gifford HC, King MA, Pretorius PH, Farncombe TH, Bruyant P, Wernick MN. Optimization of iterative reconstructions of Tc99m cardiac SPECT studies using numerical observers. IEEE Transaction on Nuclear Science. 2002;49(5):2355–2360.
  • Narayanan MV, King MA, Wernick MN, Byrne CL, Soares EJ, Pretorius PH. Improved image quality and computation reduction in 4D reconstruction of cardiac-gated SPECT images. IEEE Transactions on Medical Imaging. 2000;19(5):423–433. [PubMed]
  • Ogawa K, Harata Y, Ichihara T, Kubo A, Hashimoto S. A practical method for position-dependent Compton-scatter correction in single photon emission CT. IEEE Transactions on Medical Imaging. 1991;10:408–412. [PubMed]
  • Ogawa K, Ichihara T, Kubo A. Accurate scatter correction in single photon emission CT. Ann. Nucl. Med. Sci. 1994;7:145–150.
  • Peter J, Jaszczak RJ, Hutton BF, Hudson HM. Fully adaptive temporal regression smoothing in gated cardiac SPECT image reconstruction. IEEE Transaction on Nuclear Science. 2001;48(1):16–23.
  • Pretorius PH, King MA, Tsui BMW, LaCroix KJ, Xia W. A mathematical model of motion of the heart for use in generating source and attenuation maps for simulating emission imaging. Med. Phys. 1999;26:2323–2332. [PubMed]
  • Segars WP. PhD thesis. The University of North Carolina; 2001. Development of a new dynamic NURBS-based cardiac-torso (NCAT) phantom.
  • Zaidi H, Hasegawa B. Determination of the attenuation map in emission tomography. J. Nucl. Med. 2003;44:291–315. [PubMed]