Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Nucl Sci. Author manuscript; available in PMC 2010 April 22.
Published in final edited form as:
IEEE Trans Nucl Sci. 2009 October 1; 56(5): 2739–2749.
PMCID: PMC2858434

Theoretical and Numerical Study of MLEM and OSEM Reconstruction Algorithms for Motion Correction in Emission Tomography

Joyoni Dey, Member, IEEE and Michael A. King, Senior Member, IEEE


Patient body-motion and respiratory-motion impacts the image quality of cardiac SPECT and PET perfusion images. Several algorithms exist in the literature to correct for motion within the iterative maximum-likelihood reconstruction framework. In this work, three algorithms are derived starting with Poisson statistics to correct for patient motion. The first one is a motion compensated MLEM algorithm (MC-MLEM). The next two algorithms called MGEM-1 and MGEM-2 (short for Motion Gated OSEM, 1 and 2) use the motion states as subsets, in two different ways. Experiments were performed with NCAT phantoms (with exactly known motion) as the source and attenuation distributions. Experiments were also performed on an anthropomorphic phantom and a patient study. The SIMIND Monte Carlo simulation software was used to create SPECT projection images of the NCAT phantoms. The projection images were then modified to have Poisson noise levels equivalent to that of clinical acquisition. We investigated application of these algorithms to correction of (1) a large body-motion of 2 cm in Superior-Inferior (SI) and Anterior-Posterior (AP) directions each and (2) respiratory motion of 2 cm in SI and 0.6 cm in AP. We determined the bias with respect to the NCAT phantom activity for noiseless reconstructions as well as the bias-variance for noisy reconstructions. The MGEM-1 advanced along the bias-variance curve faster than the MC-MLEM with iterations. The MGEM-1 also lowered the noiseless bias (with respect to NCAT truth) faster with iterations, compared to the MC-MLEM algorithms, as expected with subset algorithms. For the body motion correction with two motion states, after the 9th iteration the bias was close to that of MC-MLEM at iteration 17, reducing the number of iterations by a factor of 1.89. For the respiratory motion correction with 9 motion states, based on the noiseless bias, the iteration reduction factor was approximately 7. For the MGEM-2, however, bias-plot or the bias-variance-plot saturated with iteration because of successive interpolation error. SPECT data was acquired simulating respiratory motion of 2 cm amplitude with an anthropomorphic phantom. A patient study acquired with body motion in a second rest was also acquired. The motion correction was applied to these acquisitions with the anthropomorphic phantom and the patient study, showing marked improvements of image quality with the estimated motion correction.

Index Terms: Expectation maximization algorithm, image reconstruction

I. Introduction

Patient body motion and respiratory motion for heavily breathing patients combined with the “respiratory creep” of the heart impacts the diagnostic accuracy of cardiac PET or SPECT perfusion images. In this work the focus is on motion correction, assuming the motion is known or has been estimated. A number of investigators have worked on algorithms to correct patient motions for emission tomography [1]–[12]. Some correction methods are post-reconstruction correction where different states are reconstructed individually, each of them motion corrected and then added together [3], [4], [9]. Other methods perform correction within the reconstruction algorithm [1], [2], [6]–[8], [10]. In [1], [2], [7], [8] ordered-subset expectation-maximization (OSEM) or maximum-likelihood expectation maximization (MLEM) methods have been used to correct for motion. However the EM derivations for motion correction and the assumptions made for the derivation are not explicitly shown. In [7], the proposed solution is shown to numerically increase the log-likelihood of the measurement data. Schumacher et al. [6] minimize least-squared error with a non-zero constraint to derive their formulation for motion-correction. In Schumacher’s implementation the motion is associated with the system matrix and not the object.

In this work we derive three EM motion correction algorithms from first principles assuming Poisson statistics with all the approximations and assumptions explicitly mentioned. In the first algorithm which we call Motion-Compensated MLEM (MC-MLEM), we start with the expectation of a likelihood function to maximize and after assumptions of negligible interpolation error and linear-transformation of intensities we arrived at the motion correction MLEM algorithms used in previous work [7], [8]. In doing so we also assume that the motion is invertible. For the second algorithm which we call motion-gated OSEM (MGEM-1) because we first envisioned it for use with reconstruction of respiratory-gated motion states or gates, we use the motion states as subsets in OSEM reconstruction. That is we move the measured counts (the update factor) for each of the N motion states to a fixed (reference) position to update the object at that position. Thus during each iteration we step through N subsets, updating the object N times. In a previous conference record [11], we presented another motion compensation algorithm, again using motion states as subsets. We call this method MGEM-2 herein to differentiate it from the previous method. In this method we move the object from one position to the next successively and return it to the first position only at the end of the complete iteration These three motion compensation methods are compared for bias, bias-variance characteristics and, where appropriate, the susceptibility to interpolation errors using NCAT phantom and SIMIND to simulate the SPECT data. We also applied the motion correction algorithm to a physical anthropomorphic phantom acquisition with motion simulating a steady 2 cm amplitude respiratory motion, and a clinical patient study with body motion.

II. Theory

A. Binning the Data Into Motion States

Shepp and Vardi [13] and Lange and Carson [14] have previously derived MLEM iterative reconstruction for emission tomography. In the update equation derived object motion was not explicitly taken into account. If there is motion of the object, such as patient body motion or respiratory motion, we typically have to start at an earlier stage of measurement, from the acquired list-mode data and quantize the motion into different states to form a set of projection states for the object. Within each state, the counts are assumed to be motionless. The motion between the states can be estimated by various methods, such as post-reconstruction or using the projection data themselves. The purpose of this work is to derive and evaluate different ways of reconstructing the data from these measurements. Hence we assume that the list-mode data has already been binned to different states and the motions between the states are known exactly. The starting point is therefore a set of projection data (for all the angles acquired), for each of the N motion states, Yjb, b = 0, …, N − 1. Here j is the index for j-th detector. The objective is to bring all the measured counts together to reconstruct at a reference position. In what follows, for simplicity, we index the reference state as b = 0. This reference state in actuality can be chosen to be any one of the states from 1 … N. It is understood that the indices b can be mapped in any one-to-one order to the actual states. If, for example we have 9 states, we can choose the middle State 5 as the reference state, and a typical mapping choice would be to have index states b = (0,1,2, …, 8) correspond to actual states State = (5,6,7,8,9,1,2,3,4).


Some of the common definitions and notations used to derive the three methods are explained below. In what follows M is a transformation of spatial index and Q is a motion operator operating on the intensities. This operator interpolates intensities of the object onto the new grid after transformation. M0b(i) transforms the index i to go from reference pose 0 to pose b. The transformation M is assumed to be invertible. That is M is defined as, M0b(Mb0(i))=i=Mb0(M0b(i)). It is also understood by that ib=Mb0(i) is a short hand notation indicating that the motion transformation M operates on the spatial index (x, y, z) and then the result is lexicographically ordered to ib.

The operator Q is the intensity mapping due to the motion and is defined as, (Q0bf)(i)=f(M0b(i)). Similarly, (Qb0f)(i)=f(Mb0(i)). Sometimes for notational simplicity we might write the operator as, Q0b[f(i)]. It is understood that Q0b[f(i)]=(Q0bf)(i).

We associate the unobserved complete data variables, Zijb with the object at different positions. Zijb is the count of the photons received during pose (or state) b, at detector j, that come from activity in i-th voxel (at pose b). Note that the measurement projections at each state are the sum of the random variable along the j-th projection-ray.

We define Wijb to be the random variables (RVs) obtained when we apply the motion operator on Zijb, to correspond to the activity mapped back to reference position. Thus, we have, Wijb=(Qb0Zb)ij=ZMb0(i),jb.

B. Motion Compensated MLEM Algorithm (MC-MLEM)

The expected value of the random variable Zijb, at the i-th voxel, (as seen by detector j) for the b-th state is the value of the object at the b-th transformed state, multiplied by the probability cij. In other words, cijE(Zijb)=cij(Q0bf)(i)=cijf(M0b(i))where cij is the time-interval times the probability that the i-th emission reaches the j-th detector and f(i) is the activity at the i-th location for the 0-th (reference) state. As defined in previous section, Q0b is the motion operator going from zero-th state to state b, which essentially interpolates the intensities of the object on the new grid after transformation. We formulate the expected values of the RVs Wijb=(Qb0Z)ijb corresponding to the activity transformed back to state 0. Zijb are Poisson RVs, being different noise-realizations of the object across different states and can be assumed independent across b, i, and j. The RVs Wijb are weighted sum of Poisson variables and in general not Poisson distributed. However, under the assumption that the interpolation error is small, we assume that the motion operator is mainly a rearrangement of the object in space and then would preserves the Poisson nature of the RVs in going from Z to W. In order to reconstruct the object, we wish to maximize the expection of the sum of all the RVs, across all the states, W=bWijb. The expected value of RVs Wijb are given by E(Wijb)=E((Qb0Z)ijb)cibjf(Mb0(M0b(i))=aijbf(i), where, aijb=cibj and, ib=Mb0(i).

Since we want to reconstruct the object at the reference position (position of State 0), we maximize the conditional expectation of the W, given the previous estimate of object (at reference state) and the measurements across all the states. Since W is approximately Poisson distributed, the conditional expectation we maximize is given by

E[ln(h(W,f^n+1)[mid ]Y,f^n]=bj[set membership]Sbi[set membership]lj[aijbf^n+1(i)+Lijbln(aijbf^n+1(i))+R]


Lijb=E[Wijb[mid ]Y,f^n]=E[(Qb0Z)ijb)[mid ]Yb,f^n]

using what we know of the original random variables, Z. Then we make a further approximation

LijbQb0[E(Zijb[mid ]Yb,f^n)].

This interchange between the expectation and the motion operator is possible when the motion operator is linear (weighted sum). It is not valid if the motion operator involves non-linear mapping of the variables.

But we know that E[Zijb[mid ]Yb,f^n]=cijf^n(M0b(i))Yjb/icijf^n(M0b(i)). Then it follows from (2) that


We can make one further approximation which is possible provided we ignore interpolation errors and end-effects of movement of parts of the object out of region of interest and other parts to within the region of interest. Considering the affine-transformation of two objects, if we ignore interpolation and other numerical errors, it is approximately equivalent to rotate/translate/scale/and-or/shear each of two objects by same amount, and then add (or multiply) them together, as to add (or multiply) them together first and then rotate/translate/scale/shear them. Therefore


Taking the derivative of (1) wrt fn+1, and realizing Q operates on index i and not j, and again extracting the motion operator under the assumptions mentioned above, we can further approximate, and obtain

fn+1(i)=fn(i)bQb0[j[set membership]SbcijYjbicijfn(M0b(i))]bQb0[j[set membership]Sbcij].

Notational differences apart, (4) is equivalent to the form that was used by Li et al. in [7] and Feng in [8] as solutions.

Note as expected in (4) before comparing with the measurement at b-th state, Yjb, the estimate fn (at the pose of the reference state) has to be transformed forward to the position of the b-th state, and then forward projected. Then the back-projection of the ratio of the acquired projection to the estimated projections is transformed back to 0th state and multiplied to update the previous estimate of the object at the 0th state (after normalizing with the transformed weight matrix).

C. Motion Gated OSEM Algorithm 1 (MGEM-1)

In this second approach instead of maximizing the conditional expectation of the total of the complete data variables across all the states, we maximize the expectation of the random variables at each step, conditional to the data obtained at the previous step. For notational purpose, we introduce a super-index now to correspond to states. And we introduce a sub-index (0) to under-line the fact that the object is always reconstructed to the reference position. Thus for the n-th iteration, f0b,n is the reconstructed object in the b-th state. To consider the case within a single iteration we temporarily drop the n-index.

We assume with the approximations outlined in the previous section, that the random variables Wijb are independent and Poisson. The sum of independent Poisson distributed RVs are Poisson, hence the complete data log-likelihood function is given by

ln(h(Wb,f^0b))=ji[E(Wijb)+WijblnE(Wijb)ln(Wijb!)]=ji[set membership]lj[aijbf^0b(i)+Wijbln(aijbf^0b(i))ln(Wijb!)].

We want to find the conditional expected value of the log-likelihood function of ( Wijb+1,f^0b+1) given we have the reconstruction for the previous state, f^0b, and the measurements for this state, Yb+1, and then maximize that expected value. In order to accomplish this, we first re-write the (5a) for b+1, and then take the conditional expectation, as follows:

E[ln(h(Wijb+1,f^0b+1)[mid ]Yb+1,f^b]=ji[set membership]lj[aijb+1f^0b+1(i)+Pijb+1ln(aijb+1f^0b+1(i))R]

where Pijb+1=E[Wijb+1[mid ]Yb+1,f^0b(i)] and R is the conditional expectation of the factorial which, as in [14], we need not consider further in maximization of this expectation since it does not depend on f^0b+1. But

Pijb+1=E[Wijb+1[mid ]Yb+1,f^0b]=E[(Qb+10Z)ijb+1)[mid ]Yb+1,f^0b]Qb+10(E[Zijb+1[mid ]Yb+1,f^0b]).

Again the expectation and the motion operators are interchangeable only if the motion operation is linear. Taking the derivative of (5b) wrt f^0b+1 and following similar approximations as done in the previous section, we obtain the update equation


In this equation, before comparing with new measurements at the (b + 1) state, Yjb+1, the previous estimate f^0b has to be transformed to the position of the (b + 1) state and then forward projected. The back-projection of the ratio is transformed back to 0th state and multiplied to update the object (after normalizing with the transformed weight matrix). After stepping through all the states, we can iterate the process.

D. Motion Gated OSEM Algorithm 2 (MGEM-2)

The second version of the MGEM algorithm differs from the first by how the motion states used as subsets are employed to update the estimate of the source distribution. Here we directly work with the random variables Zijb. For notational clarity, we assume that the different states indexed by b provide information on the transformed versions of the object, f, given by transformed states, f0, …, fb, …, fN−1. This is consistent with previous section’s notation, where we dealt with f0(i), the object at the reference position. The transformed states are defined by fb=Q0b[f0(i)].

The expected value of the random variable Zijb, is given by


We then reconstruct the object f using intermediate steps of reconstructing the transformed states f0, …, fb, …, fN − 1 and use one final transformation to go from the last amplitude state to the first before iterating the process. Let Tb+1 be the transformation operator to take the b-th state to the (b + 1)th state. This Tb+1(fb) is a transformation of the indices i and interpolation of the intensities on the transformed grid. TN is the transformation to go from the N − 1 back to the 0-th state.

We want to find the conditional expected value of the log-likelihood function of (Zb+1, fb+1) given we have the previous reconstructed stage, fb, and the new measurements, Yb+1, and then maximize that expected value. Hence,

E[ln(h(Zb+1,f^b+1)[mid ]Yb+1,f^b]=ji[set membership]lj[cijf^b+1(i)+Dijb+1ln(cijf^b+1(i))R]

where Dijb+1=E[Zijb+1[mid ]Yb+1,f^b(i)] and R is the conditional expectation of the factorial which, as in [14], we need not consider further in maximization of this expectation since it does not depend on fb+1. But,

Dijb+1=E[Zijb+1[mid ]Yb+1,f^b]=cijTb+1(f^b)Yjb+1i[set membership]ljcijTb+1(f^b).

This is explained in two steps. First, our previous estimate fb is at a different location and therefore needs to be transformed to give the estimate of the volume at the current location. Hence, the fb is transformed to the next state’s position by Tb+1(fb) first. Then, as in [14], we can derive the expected value at the i-th pixel (and j-th detector) by dividing the new measurement Yjb+1 into the ratio of the voxel estimate that that location by the summed expected values along the ray. Thus we get our Dijb+1 as shown above. Taking the derivative of (6), with respect to fb+1, and rearranging, we would get the update equation as,


where again, Tb+1 is the transformation from position of state b to b + 1.

After one pass through all the states, we make sure to transform to the first state.


We can iterate this process and introduce an iteration variable k as a superscript.

Our general update Eqns will be given by


At the end of this iteration, after going through all the states b = 0 …, N − 1, we go back to the zero-th state as


Then we go on to the next iteration as


Equations (8)(9)(9a) together shows our update scheme. A clear advantage of (8) over the update in (4) is that in using (8), only one transformation of the object is necessary for each stage while for (4) to obtain the update factor, we need to forward transform the object to current state from reference state (before forward projection to compare to measurement) and then after back-project apply an inverse transform to bring it back to reference. Also the linearity of intensity transforms is not assumed explicitly. However a major drawback of (8) is in reality the object would be interpolated with each stage and also for each iteration. This could affect the smoothness of the reconstruction.

In what follows we perform experiments on NCAT phantom, an anthropomorphic phantom and a patient study to compare these motion compensation methods.

III. Methods

A. Study of Rigid-Body-Motion Correction Using NCAT Simulated Acquisitions

To simulate body-motion NCAT phantoms were generated for the source and attenuation distribution in two different motion states with the first displaced from the second with 20 mm of motion in Z direction and 20 mm motion in Y direction. For each of these NCAT datasets, SIMIND Monte-Carlo software was used to simulate clinical SPECT acquisition for Tc-99m sestamibi for 60 angular projections acquired over 180 degrees at steps of 3 degrees. Solely the primary events simulated by SI-MIND are employed herein. Thus scatter is not included in the photopeak window and scatter correction is not included in reconstruction. It was assumed that the acquisition was acquired by a 2-headed SPECT system, thus of the 60 angular projections, first 30 belonged to one head and next 30 belonged to the second-head. Then these projections from the two datasets were used selectively to simulate acquisition for the patient in two rigid-body motion states, with the projection sets, Yj0 and Yj1. The projections for the two motion states were then merged to create a projection set in which there is no motion in the first half of the acquisition (angles 1 to 15 and angles 31 to 45) and then 2 cm motion (in Y and Z) for the second half of the acquisition (angles 16 to 30 and 46 to 60). Random Poisson noise was added to the projections such that the total count over two states was 7.5 million. 50 noise samples of the projection states were generated and the object was reconstructed with the MC-MLEM, MGEM-1 and MGEM-2 for each of the 50 cases, as well as for the noiseless projections. Reconstruction included attenuation correction and modeling of the distance-dependent system spatial-resolution [15]. The attenuation map for the first motion state was employed in reconstruction and it was moved using the known motion between the two states prior to projection and backprojection for projections acquired from the second motion state.

For the purpose of evaluation of our motion compensation algorithms, we observe how closely the outputs of the motion compensating algorithms compare with the motionless object, that is, the original NCAT phantom activity distribution for the first motion state (reference state). Thus the reconstructions from projections without added Poisson noise were compared voxel-wise with the NCAT at the reference stage. The bias was calculated as the absolute difference of the activity of the NCAT phantom from the mean of the reconstructed object over a region of the heart, normalized by the average of the reference over the region (NCAT). We also investigated the bias-variance characteristics for reconstructions with the reconstructions of the 50 noise realizations. The bias was calculated as the absolute difference of the activity of the NCAT phantom from the ensemble mean (over 50 noise realizations), averaged over a region of the heart, and normalized with respect to the NCAT mean. The variance was calculated with respect to the ensemble mean and then averaged over the region of the heart. Bias-variance plots were calculated as a function of iteration number to compare the no motion case, the motion-case without correction and the motion correction cases with our 3 reconstruction methods. We also looked into tri-linear interpolation versus Gaussian interpolation for the different algorithms.

B. Study of Respiratory Motion Correction Using NCAT Simulated Acquisitions

Thirty-six NCAT phantoms were generated to simulate a respiratory-amplitude binned acquisition for a Tc-99m sestamibi imaging study. In these source distributions and attenuation maps the respiratory motion of the heart equally dividing 20 mm in Z and 6 mm in Y net motion of the heart. SIMIND Monte-Carlo projections of each were generated to simulate clinical SPECT acquisition. Each of 4 projection sets were added at a time to create 9 states out of the 36 states. This was to make sure we have realistic effects of finite motion-quantization. That is, each of 9 bins would not have exactly stationary counts, but have an average of 4 motion states. Poisson noise was added such that total count (over all 9 bins) over the torso was 7.5 million. Fifty such noise-instances were considered. Each sample of the object was then reconstructed using the MC-MLEM, MGEM-1 and MGEM-2 algorithms, choosing the middle-state as the reference state. Reconstruction included attenuation correction and modeling system spatial resolution as before.

Considering (4), we observe that the update factor has a summation over the states b in the numerator and the denominator. For respiratory motion considered here, the states are comprised of all the angles. For the noisy cases, the projection set for each state would have 1/9 the number of the counts of the total (ie, 7.5/9 = 0.833 million each). In other words, the counts in each angle at each state will have a fraction of counts acquired for that angle. In (4), the summation over states b in the numerator would yield an overall sum of all counts in all the states, but there would be an effective factor 9 in the denominator as well, from the addition of the 9 sensitivity matrices. Thus the counts in the resultant reconstruction would be expected to be scaled down by a factor of 9 (in general, by the no. of states). Similarly, since the counts at each angle are a fraction of total acquired at the angle, the MGEM algorithms for the noisy datasets for respiratory motion will also be driven to a net count of 1/9 of the total counts. Hence, for the noise-added cases, the objects reconstructed with the motion correction respiratory motion algorithms were scaled appropriately before bias-variance comparison to no motion case or the no correction case (where there was motion but no correction was applied).

Note that for body motion the motion states are comprised of mutually exclusive angles and each angle has the full acquired number of counts. Thus the reconstructed counts are preserved to be the total counts in all the states.and this scale factor was not an issue to consider for body motion.

As for the body-motion case, the bias was calculated from the noiseless reconstructions with the NCAT phantom. The reference state was the middle one of the 9. Bias-variance plots were created from analysis of the 50 noisy reconstructions for a region around the heart. Due to the scale factor issue as explained above, we scaled each of the datasets by the ensemble means averaged over the region considered, for a fair comparison between the algorithms. On numerical comparison, the ratio of the scale factors for the no-motion case and the motion correction cases was indeed found to be close to 9. This scale factor issue was taken into account for the bias part for the bias-variance calculation, as well – the NCAT and the data were each scaled by their respective averages over region of the heart before taking the absolute difference between them.

C. Study of Respiratory Motion Correction Using Acquisitions of the Data Spectrum Anthropomorphic Phantom

The Data Spectrum Anthropomorphic phantom with Tc-99m in the heart wall, liver, and background was imaged with an IRIX SPECT system. List-mode data was acquired with the phantom moved manually in a periodic fashion to simulate respiratory motion with an amplitude of 2 cm. The motion was regular, in approximate synchrony with the breathing of the person moving the phantom. Motion was monitored using a bellow attached to the phantom and the imaging table [11]. A second study was acquired where the data was not moved for purposes of comparison. For the motion-case the acquired list-mode data was binned into 9 different amplitude states. To estimate the motion, each of the states was then reconstructed and the motion was estimated using an intensity-based registration method. The registration uses sum-squared-error as the similarity measure to minimize and iteratively uses gradient descent for the minimization [16]. The method allows for 1 to 12 degrees of freedom (DOF) motion estimation. For the case in consideration here, even though experimentally there was 1-DOF motion, 2-DOF – Superior-Inferior (SI) and anterior-posterior (AP) – estimates were performed as these are the primary directions of expected motion for the heart with respiration [16], [17]. The motion estimates and the amplitude-binned projection states were used for the reconstruction correction methods described herein. The motion corrected datasets were compared to the case where there was no motion correction case as well as the no-motion case. Again reconstruction included attenuation correction and modeling of distance-dependent spatial resolution. The attenuation map was acquired at the no-motion case (without simulated respiratory motion) and was used in correction of the studies which included respiratory motion, by appropriately moving the map using the motion estimates.

D. Rigid-Body Motion Correction for a Patient Study

We have developed a stereo-imaging based motion-tracking system which estimates the motion of the heart from the motion of retro-reflective markers on bands wrapped about the chest and abdomen of patients [18]. With Institutional Review Board approval and informed consent motion-tracking was performed in a series of patients undergoing two T1-201 cardiac-perfusion SPECT rest-imaging acquisitions on our IRIX SPECT system. During the first acquisition the patients remain as motionless as possible as standard during SPECT imaging. During the second acquisition the patients performed one or more motions at our request. This made available the initial rest-imaging study to provide a reference to compare the appearance of the second study with successful motion estimation and correction. For purpose of illustration of our methods herein, we chose a patient who slid down the table about ~ 6 cm axially. The translation of this patient was determined in the axial, lateral, and vertical directions from the motion tracking-data. We then used the estimated motion and the acquired projections to reconstruct using the MC-MLEM and MGEM-1 methods described here. The motion corrected datasets were compared with the first rest study. For the MGEM-1 the data was divided into two motion angular subsets as there was a clear demarcation in our experiment. Reconstruction included attenuation correction using transmission maps reconstructed by transmission profiles acquired using the Beacon system between two rest-imaging studies. Again, the map was moved using the motion estimates when used for attenuation correction for the motion correction algorithms.

IV. Results

A. Study of Rigid-Body-Motion Correction Using NCAT Simulated Acquisitions

In Fig. 1 noise-free short-axis cardiac SPECT slices are shown before and after motion correction. In all cases reconstruction included attenuation correction and correction for distance-dependent spatial resolution. Scatter correction is not included as scatter was not included in the projection sets. Fig. 1(a) shows the case where no motion correction was applied when motion was simulated as being present. Standard MLEM reconstruction [13], [14] was employed and 18th iteration is shown. Notice the significant degradation of image quality caused by not correcting for motion. This degradation can be qualitatively appreciated by comparison to Fig. 1(b) which shows the results of 18 iterations of MLEM reconstruction for simulated projections form the phantom just in the reference state. Fig. 1(c) shows a motion corrected case using MC-MLEM after 18 iterations, which qualitatively looks similar to the reconstruction for the phantom just in the reference state. Fig. 1(d)–(e) show the reconstruction of the projections included motion using MGEM-1 after 9 iterations and 18 iterations respectively. Qualitative comparison seems to indicate that the resolution at iteration 9 of MGEM-1 is close to that of iteration 18 of MC-MLEM, and 18 iterations of MGEM-1 shows even more resolution recovery. This is as expected since there are two motion states, with mutually exclusive angles, thus this is an OSEM reconstruction with 2 angular subsets. Fig. 1(f) shows the results of 18 iterations of MGEM-2 reconstruction. It shows an unexpected amount of smoothing. For example the 18th iteration of MGEM-2 shows more smoothing than the 18th iteration of MC-MLEM in Fig. 1(c) even though MGEM-2 should also act like an OSEM reconstruction with 2 subsets. To investigate the reason for this, we created another NCAT dataset with motion exactly 4 pixel of motion (or 1.87 cm) in Y and Z as opposed to the 2 cm motion in the previous simulations. This ensured that there would be no interpolation error as the motion is an integer number of pixels. With this interpolation-error-free case, the behavior of MGEM-2 was similar to MGEM-1, as shown in Fig. 1(g). These results Fig. 1(f)–(g) indicates that (8) is prone to successive interpolation error and therefore the MGEM-2 method should be avoided.

Fig. 1
NCAT Body-Motion Correction of noiseless datasets. (a) Mid-ventricular short-axis slice from simulation which did include motion with no motion correction included in 18 iterations of MLEM. This slice illustrates the significant degradation uncorrected ...

Fig. 2 shows the normalized bias for the noise-free reconstructions with respect to the reference NCAT phantom for 100 iterations. The MGEM-1 reconstruction was observed to be lowering its bias with respect to the NCAT phantom at a faster rate than the MC-MLEM, as expected. MGEM-1 being a subset-type algorithm is expected to reduce the number of iterations by a factor of number of subsets (the number of motion bins in this case). The bias value for the MGEM-1 case at iteration ~ 9 was close to that of the MLEM (tri-linear) case at iteration 17, indicating an iteration reduction factor of 1.89 around this operating point, which is close to the number of states (2). For MC-MLEM, the tri-linear interpolation seems to have a lower bias than the Gaussian interpolation. The bias values of MGEM-2 (tri-linear) on the other hand seem to flatten, qualitatively speaking. Bias of the MGEM-2 Gaussian actually increases before flattening. However when we consider the bias for the integer motion case (MGEM-2-int), the behavior is close to MGEM-1, indicating that, indeed, the higher error was due to interpolation. Another interesting point is that the no-motion case shows bias similar to as MC-MLEM for the first 20 iterations and then gradually converges faster until it is more like an MGEM algorithm.

Fig. 2
NCAT Body-Motion Correction: Normalized Bias as a function of iteration for various reconstruction methods for noiseless NCAT projection sets. Note that (G) signifies the use of 3D Gaussian interpolation, (T) signifies the use of trilinear interpolation, ...

Fig. 3 shows a bias-variance plot for the various reconstruction algorithms for 100 iterations in each case. Again the MGEM-1 advanced along the bias-variance plot faster than MLEM. Also, MGEM-2 with non-integer motion did not behave like an OSEM algorithm except with integer motion (MGEM-2-int). In fact in presence of noise the MGEM-2-int (motion correction without interpolation error) had initially similar performance as MGEM-1 but gradually lowered the bias more with iterations. Note that the no-motion MLEM case lowered the bias furthest.

Fig. 3
NCAT Body-Motion Correction: Normalized-Bias vs. Variance of Ensemble Mean of 50 noisy reconstructions for various reconstruction methods for noiseless NCAT projection sets. Note that (G) signifies the use of 3D Gaussian interpolation, (T) signifies the ...

B. Study of Respiratory Motion Correction Using NCAT Simulated Acquisitions

Fig. 4(a) shows an NCAT short-axis slice for 18 iterations of MLEM reconstruction of noise-free projections containing respiratory-motion without motion correction. The significant degradation of the respiratory motion can be appreciated by comparing this slice to reconstruction of projections from just the reference state (5th state) in Fig. 4(b). Notice the increase apparent ventricular wall thickness, decrease in the size of the blood-pool region, and increased blurring with the liver in the presence of respiratory motion. The motion correction algorithms as shown in Fig. 4(c)–(e) clear up these degradations. Qualitative comparisons between the slices of Fig. 4(c) and (d) seem to indicate that the resolution at iteration 2 of MGEM-1 is close to that of iteration 18 of MC-MLEM, as expected here with 9 respiratory states. In Fig. 4(f) we show MGEM-2. Again there is uncharacteristic smoothing assumed to be due to interpolation error. Since 1 pixel integer motion corresponding to each bin will result in a large motion, we did not test the interpolation-error-free case for this respiratory motion.

Fig. 4
NCAT Respiratory motion correction of noiseless datasets. (a) Mid-ventricular short-axis slice from simulation which did include respiratory motion with no motion correction included in 18 iterations of MLEM. (b) Short-axis slice from projections simulated ...

Fig. 5 shows the normalized bias plots for the noise-free reconstructions for up to 100 iterations of the reconstruction algorithms, using the NCAT activity as a comparison. Again we observe that the MGEM-1 seem to achieve closer values to the NCAT faster than the MLEM algorithms. In particular, the bias value at iteration 2 of MGEM-1 is similar to the bias value at iteration 14 for MC-MLEM (with tri-linear interpolation), indicating a factor of 7 reduction in iterations around this operating point for 9 subsets in MGEM-1. Again the no-motion case performs close to MC-MLEM in the first ~ 35 iterations and then gradually lowers the bias faster. The MGEM-2 with Gaussian (or Trilinear) interpolation seems to lower the bias for only a few iterations and then levels it out.

Fig. 5
NCAT Respiratory motion Correction: Normalized Bias as a function of iteration for various reconstruction methods for noiseless NCAT projection sets. Note that (G) signifies the use of 3D Gaussian interpolation and (T) signifies the use of trilinear interpolation. ...

Fig. 6(a) shows the bias-variance plot for 100 iterations, derived as explained in the methods section. MGEM-1 (with 9 motion states used as subsets) advanced significantly faster along the bias-variance curve than the MC-MLEM or the other MLEM algorithms (such as the no-motion case or the no-motion-correction cases). This was as expected due to the OSEM nature of MGEM-1, and high number of motion states. For better visualization we displayed a zoomed-in view in Fig. 6(b), which shows that the MGEM-1 achieves the similar bias-variance in first few iterations (~ 13) as the No-Motion case (MLEM algorithm) did after 100 iterations.

Fig. 6
(a) NCAT Respiratory Motion Correction: Normalized-Bias vs. Variance of Ensemble Mean of 50 noisy reconstructions for various reconstruction methods for noiseless NCAT projection sets, shown for 100 iterations points. Note that (G) signifies the use of ...

For a numerical comparison, after 100 iterations, the No-Motion case had bias variance of (0.41,0.04); MC-MLEM (G) was at (0.49,0.03); while the MGEM-1 had those at (0.42,0.47). Note the order of magnitude higher variance for the MGEM-1 (expected of a subset algorithm, proceeding faster than the MLEM ones), However, a slight increasing trend of bias was observed for later iterations of MC-MLEM (T or G) and MGEM-1. The performance of the MGEM-2, (Gaussian or Trilinear) was similar to as observed for the noiseless case, in that was slower than even the MLEM algorithms and quickly saturated.

C. Study of Respiratory Motion Correction Using Acquisitions of the Data Spectrum Anthropomorphic Phantom

The registration estimated motion across the 9 amplitude bins is −10.35, −6.92, −4.60, −1.75, 0, 3.33, 5.30, 8.7, and 9.93 mm in the superior-inferior (SI) direction which matches well the intended motion of +/− 10 mm in this direction, and 0.62, 0.36, 0.07, 0.04, 0, −0.01, 0.01, 0.07, and 0.90 mm in anterior-posterior (AP) direction which matches the intended no motion in this direction. Fig. 7 shows a mid-ventricular short-axis slice for reconstruction which in each case included attenuation and resolution compensation. For the motion corrected cases, the attenuation map obtained at the no-motion reference state was moved according to the estimated motion for the different states. Fig. 7(a) shows the 26 MLEM iterations for the case where the phantom was kept stationary on the table. This is the no-motion reference state. Fig. 7(b) shows 36 iterations of MLEM without motion correction for the case with motion (the acquisition performed with respiratory-motion created by having a volunteer move the phantom in approximate synchrony with their breathing). Notice the false cooling in the superior-inferior walls. Fig. 7(c)–(d) shows motion correction with MC-MLEM (with Gaussian interpolation) and MGEM-1, respectively. The no-motion was seen to advance faster than the MC-MLEM, and we chose to compare iteration 26 of MLEM for the no motion to 36 of MC-MLEM of the motion correction as they had similar wall thickness. MGEM-1 achieves similar results in just 4 OSEM iterations. This is roughly 9 times faster than the MC-MLEM which is expected with the 9 motion sub-sets.

Fig. 7
Data Spectrum Anthropomorphic Phantom: (a) Mid-ventricular short-axis slice after 26 iterations of MLEM for acquisition during which the phantom was not moved. (b) Short-axis slice after 36 iterations of MLEM without motion correction for acquisition ...

D. Rigid-Body Motion Correction for a Patient Study

The three translation estimates for the patient with angle of acquisition are shown in Fig. 8. Fig. 9 shows short-axis slices for the reconstructions which in each case included attenuation and resolution compensation. For the motion corrected cases, the attenuation map obtained at the initial rest reference state was moved according to the estimated motion for the different states. Fig. 9(a) shows the initial rest where the patient was instructed not to move. Fig. 9(b) shows the reconstructed second rest without motion compensation during which the patient was requested to perform an “axial-slide” after a time interval. Fig. 9(c)–(d) shows the results for MC-MLEM (Gaussian interpolation) and MGEM-1 with motion correction. The first rest was again seen to converge slightly faster in terms of wall thickness than MC-MLEM, hence we compared iterations 30 of MLEM to 40 with MC-MLEM. Again MGEM-1 achieves similar results in 20 iterations, as expected for 2-subset motion correction.

Fig. 8
Estimated translational motion of the patient as a function of acquisition angle for acquisition by two camera-heads at 102 degrees.
Fig. 9
Patient Body Motion Correction: (a) Mid-ventricular short-axis slice from initial-rest study acquired without intentional patient motion and reconstructed with 30 iterations of MLEM. (b) Short-axis slice from second-rest study acquired with intentional ...

V. Discussion

The primary aim of this work was to derive a theoretical framework for motion-correction within an MLEM iterative algorithm, (MC-MLEM) with all the approximations mentioned explicitly. The two essential theoretical findings were that the intensity transformations involved in the motion have to be linear and invertible. Also couple of fast methods of motion compensation (MGEM-1 and MGEM-2) were also derived and compared. We used the NCAT phantom as well as physical phantom and patient acquisitions to study the motion correction methods. We compared the performance of the MLEM and two MGEM (“motion-gated” OSEM, where the motion states serve as subsets) for body motion and respiratory motion. One MGEM-1 moved back the measurement counts for each motion state to the reference state position, and then the estimate at the reference state location was updated. We found expectedly that this MGEM algorithm converged faster towards the true NCAT object by a factor close to the number of states, compared to the MLEM correction method. In MGEM-2 we successively moved from one motion state to the next and then finally back to reference. While this algorithm seemed good in theory, it had successive interpolation artifacts which resulted in significant smoothing and was therefore not useful.

OSEM reconstruction is typically applied with each subset containing an equal number of projections [19]. In the past Feng, et al. [8] showed that one could include motion into OSEM by keeping the same subset structure as one would normally use by using 3D Gaussian interpolation to account for motion at any angle where needed. Herein we expanded use of motion correction by considering subsets formed by respiratory binning or guided by body motion. For respiratory motion considered here, each subset contained a portion of the counts acquired at all projection angles. For body-motion subsets were formed as being acquired at the same (body) motion state so that there may not in general be the same number of angles in each subset. Kyme, [20] have also previously formed body-motion subsets with an unequal number of projection angles in the subsets. For the patient-study we considered the results of use of an unequal number of angles was encouraging; however, as the number of angles in one or more subsets becomes very different from the others it may be that using the rescaled-block-iterative (RBI) algorithm [21] may be needed to avoid artifacts. Also as explained in the methods section, in the case where the motion correction algorithms were applied to the amplitude-binned respiratory motion with the motion-states made up of projections containing only a portion of the counts acquired at each angle, it is important to include a scale factor before comparing to no-motion case.


This work was supported by the National Institute of Biomedical Imaging and Bioengineering (NIBIB) grant R01 EB001457. Part of this work was presented at the 2008 IEEE NSS-MIC [22].

The authors would like to thank Prof. C. Byrne from University of Massachusetts, Lowell, for his insightful discussions and help with the notation of this paper. They would also like to thank Dr. J. Mitra at University of Massachusetts Medical School, Worcester, for his helpful review of this paper. The contents are solely the responsibility of the authors and do not necessarily represent the official views of the National Institutes of Health.


1. Fulton RR, Hutton BF, Braun M, Ardekani B, Larkin RS. Use of 3-D reconstruction to correct for patient motion in SPECT. Phys Med Biol. 1994;39:563–574. [PubMed]
2. Hutton BF, Kyme AZ, Lau YH, Skerrett DW, Fulton RR. A hybrid 3-D reconstruction/registration algorithm for correction of head motion in emission tomography. IEEE Trans Nucl Sci. 2002 Feb;49(1):188–194.
3. Klein GJ, Reutter RW, Huesman RH. Four-dimensional affine registration models for respiratory-gated PET. IEEE Trans Nucl Sci. 2001 Jun;48(3):756–760.
4. Livieratos L, Stegger L, Bloomfield PM, Schafers K, Bailey DL, Camici PG. Rigid-body transformation of list-mode projection data for respiratory motion correction in cardiac PET. Phys Med Biol. 2005;50(14):3313–3322. [PubMed]
5. Kovalski G, Israel O, Keidar Z, Frenkel A, Sachs J, Azhari H. Correction of heart motion due to respiration in clinical myocardial perfusion SPECT scans using respiratory gating. J Nucl Med. 2007;48(4):630–636. [PubMed]
6. Schumacher H, Fischer B. A new and flexible reconstruction framework for motion correction in SPECT Imaging. IEEE Trans Nucl Sci. 2007 Jun;54(3):480–485.
7. Li T, Thorndyke B, Schreibmann E, Yang Y, Xing L. A model-based image reconstruction for four-dimensional PET. Med Phys. 2006;33(5):1288–1298. [PubMed]
8. Feng B, Gifford HC, Beach RD, Boening G, Gennert MA, King MA. Use of three-dimensional Gaussian interpolation in the projector/backprojector pair of iterative reconstruction for compensation of known rigid-body motion in SPECT. IEEE Trans Med Imag. 2006 Jul;25:838–844. [PMC free article] [PubMed]
9. Lamare F, Carbayo MJL, Cresson T, Kontaxakis G, Santos A, Rest CCL, Reader AJ, Visvikis D. List-mode- based reconstruction for respiratory motion correction in PET using non-rigid body transformations. Phys Med Biol. 2007;52(17):5178–520. [PubMed]
10. Dey J, Feng B, Johnson KL, McNamara JE, Pretorius PH, King MA. Respiratory motion correction in cardiac SPECT using affine and free-form registration with temporal and spatial constraints. presented at the 9th Int. Meet. Fully 3D; Lindau, Germany. Jul. 2007.
11. Dey J, Johnson K, Pretorius PH, Mitra J, King MA. Impact of variability in respiration on estimation and correction of respiratory motion of the heart in SPECT. Conf. Rec., 2007 IEEE NSS-MIC; Hawaii. Oct./Nov. 2007; pp. 4503–4509.
12. Dey J, Segars WP, Pretorius PH, King MA. Estimation and correction of irregular respiratory motion of the heart in SPECT in presence of partial angle effects due to amplitude binning in SPECT. presented at the IEEE NSS-MIC; Dresden, Germany. Oct. 2008.
13. Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imag. 1982 Oct;1(2):113–122. [PubMed]
14. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assisted Tomography. 1984 Apr;:306–316. [PubMed]
15. Narayanan MV, King MA, Pretorius PH, Dahlberg ST, Spencer F, Simon E, Ewald E, Healy E, MacNaught K, Leppo JA. Human observer ROC evaluation of attenuation, scatter and resolution compensation strategies for Tc-99m myocardial perfusion imaging. J Nucl Med. 2003 Nov;44(11):1725–1734. [PubMed]
16. Dey J, Pan T, Choi DJ, Robotis D, Smyczynski M, Pretorius PH, King MA. Semi-automated segmentation of non-contrast-enhanced 4D-CT cardiac datasets with application to respiratory-motion estimation of the heart. IEEE Trans Nucl Sci. 2009 Jan; submitted for publication. [PMC free article] [PubMed]
17. McLeish K, Hill DLG, Atkinson D, Blackall JM, Razavi R. A study of the motion and deformation of the heart due to respiration. IEEE Trans Med Imag. 2002 Sep;21(9):1142–1150. [PubMed]
18. McNamara JE, Pretorius PH, Johnson K, Mitra J, Dey J, Gennert MA, King MA. A flexible multi-camera visual-tracking system for detecting and correcting motion-induced artifacts in cardiac SPECT slices. Med Phys. 2009 May;30(5):1913–1923. [PubMed]
19. Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imag. 1994 Dec;13(4):601–609. [PubMed]
20. Kyme AZ, Hutton BF, Hatton RL, Skerrett DW, Barden LR. Practical aspects of a data-driven motion correction approach for brain SPECT. IEEE Trans Med Imag. 2003 Jun;22(6):722–729. [PubMed]
21. Byrne CL. Accelerating the EMML algorithm and related iterative algorithms by rescaled block-iterative methods. IEEE Trans Imag Process. 1998 Jan;7(1):100–109. [PubMed]
22. Dey J, King MA. Theoretical and numerical study of MLEM and OSEM reconstruction algorithms for motion correction in emission tomography. presented at the Conf. Rec. IEEE NSS-MIC; Dresden, Germany. Oct. 2008. [PMC free article] [PubMed]