|Home | About | Journals | Submit | Contact Us | Français|
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.
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 –. Some correction methods are post-reconstruction correction where different states are reconstructed individually, each of them motion corrected and then added together , , . Other methods perform correction within the reconstruction algorithm , , –, . In , , ,  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 , the proposed solution is shown to numerically increase the log-likelihood of the measurement data. Schumacher et al.  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 , . 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 , 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.
Shepp and Vardi  and Lange and Carson  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, , 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. 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, . It is also understood by that 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 .
The operator Q is the intensity mapping due to the motion and is defined as, . Similarly, . Sometimes for notational simplicity we might write the operator as, . It is understood that .
We associate the unobserved complete data variables, with the object at different positions. 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 to be the random variables (RVs) obtained when we apply the motion operator on , to correspond to the activity mapped back to reference position. Thus, we have, .
The expected value of the random variable , 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, 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, 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 corresponding to the activity transformed back to state 0. 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 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, . The expected value of RVs are given by , where, and, .
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
using what we know of the original random variables, Z. Then we make a further approximation
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 . 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 n+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
Note as expected in (4) before comparing with the measurement at b-th state, , 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).
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, 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 are independent and Poisson. The sum of independent Poisson distributed RVs are Poisson, hence the complete data log-likelihood function is given by
We want to find the conditional expected value of the log-likelihood function of ( ) given we have the reconstruction for the previous state, , 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:
where and R is the conditional expectation of the factorial which, as in , we need not consider further in maximization of this expectation since it does not depend on . But
Again the expectation and the motion operators are interchangeable only if the motion operation is linear. Taking the derivative of (5b) wrt 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, , the previous estimate 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.
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 . 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 .
The expected value of the random variable , is given by
We then reconstruct the object f using intermediate steps of reconstructing the transformed states 0, …, b, …, N − 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(b) 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, b+1) given we have the previous reconstructed stage, b, and the new measurements, Yb+1, and then maximize that expected value. Hence,
where and R is the conditional expectation of the factorial which, as in , we need not consider further in maximization of this expectation since it does not depend on b+1. But,
This is explained in two steps. First, our previous estimate b is at a different location and therefore needs to be transformed to give the estimate of the volume at the current location. Hence, the b is transformed to the next state’s position by Tb+1(b) first. Then, as in , we can derive the expected value at the i-th pixel (and j-th detector) by dividing the new measurement into the ratio of the voxel estimate that that location by the summed expected values along the ray. Thus we get our as shown above. Taking the derivative of (6), with respect to b+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.
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, and . 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 . 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.
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.
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 . 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 . 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 , . 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.
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 . 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.
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 ,  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. 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. 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. 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. 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. 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.
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.
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.
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.
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 . In the past Feng, et al.  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,  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  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 .
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.