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 2010 July 14.
Published in final edited form as:
PMCID: PMC2903972

LOR-OSEM: statistical PET reconstruction from raw line-of-response histograms


Iterative statistical reconstruction methods are becoming the standard in positron emission tomography (PET). Conventional maximum-likelihood expectation-maximization (MLEM) and ordered-subsets (OSEM) algorithms act on data which has been pre-processed into corrected, evenly-spaced histograms; however, such pre-processing corrupts the Poisson statistics. Recent advances have incorporated attenuation, scatter, and randoms compensation into the iterative reconstruction. The objective of this work was to incorporate the remaining preprocessing steps, including arc correction, to reconstruct directly from raw unevenly-spaced line-of-response (LOR) histograms. This exactly preserves Poisson statistics and full spatial information in a manner closely related to listmode ML, making full use of the ML statistical model. The LOR-OSEM algorithm was implemented using a rotation-based projector which maps directly to the unevenly-spaced LOR grid. Simulation and phantom experiments were performed to characterize resolution, contrast, and noise properties for 2D PET. LOR-OSEM provided a beneficial noise-resolution tradeoff, outperforming AW-OSEM by about the same margin that AW-OSEM outperformed pre-corrected OSEM. The relationship between LOR-ML and listmode ML algorithms was explored, and implementation differences are discussed. LOR-OSEM is a viable alternative to AW-OSEM for histogram-based reconstruction with improved spatial resolution and noise properties.


Iterative statistical reconstruction methods are becoming the standard in positron emission tomography (PET) and single-photon emission computed tomography (SPECT). These methods model the statistical nature of the measured data, e.g. Poisson statistics, and they can be classified into two broad categories: histogram-based methods, which reconstruct from binned histograms of the measured data; and listmode methods, which reconstruct from event-by-event lists of measured events. Statistical reconstruction of histogrammed PET data is largely based on the maximum-likelihood expectation-maximization (MLEM) algorithm described in the papers by Shepp and Vardi (1982) and Lange and Carson (1984). The MLEM update equation can be written:


where x̂i is the value of voxel i of the reconstructed image (i=1…N), [p with tilde]j is the (noisy) number of prompt coincidence events measured in projection histogram bin j (j=1…M), Fji is the probability that an emission within voxel i gives rise to a coincidence measurement at projection bin j, and n is the iteration number.

In its purest form, each histogram bin of raw PET data exactly corresponds to a single line-of-response (LOR) between a pair of detector elements. For a ring tomograph, the LORs are unevenly spaced in s as shown in figure 1. Most conventional implementations of MLEM and related algorithms interpolate the raw LOR histogram, in a process called arc correction, into a projection image matrix which is evenly spaced in s (Defrise and Kinahan 1998). This involves undesirable interpolation, but has been the historical standard largely due to the ease of reconstructing evenly-spaced data. Additional pre-processing steps are also commonly performed as shown in figure 2. When any such pre-processing is applied, it spoils the Poisson nature of the data and the ML statistical model is under-utilized.

Figure 1
Diagram of a generic ring PET tomograph showing a set of LORs (left) and the (s, [var phi], z, δ) coordinate system used to parameterize each LOR (center, right). Due to the cylindrical arrangement of crystals the LORs grouped into a parallel ...
Figure 2
Data processing schemes for conventional MLEM (top), AW-MLEM (middle), and the proposed LOR-MLEM (bottom). Note the progression from pre-correcting for image measurement effects to incorporation of these effects into the system transfer matrix used for ...

Significant developments in the field include accelerating convergence (e.g. ordered-subsets and row-action algorithms—Hudson and Larkin 1994, Byrne 1996, Browne and De Pierro 1996), regularizing the solution to reduce noise artifacts (e.g. Bayesian or penalized-likelihood algorithms), and correcting or compensating for image degrading factors (attenuation, scatter, randoms, point response). A recent improvement particularly pertinent to this work is the so-called attenuation-weighted ordered-subsets (AW-OSEM) algorithm (Comtat et al. 1998, Hebert and Leahy 1990). In AW-OSEM, the process of non-uniform attenuation is explicitly modeled in F as opposed to pre-correcting the data. This has the advantage of keeping the data “more Poisson-like”, and it results in a more statistically efficient reconstruction which produces lower noise images. Indeed, Shepp and Vardi (1982) recognized that all measurement effects should be modeled in F in order to make full use of the ML statistical model.

Unlike the histogram-based methods, listmode methods store an event-by-event list of the raw PET data attributes (e.g. position and energy) as the data are acquired. The listmode ML algorithm reconstructs directly from the list without binning into a histogram (Barrett et al. 1997, Parra and Barrett 1998, Reader et al. 1998, Huesman et al. 2000, Bouwens et al. 2001, Byrne 2001, Levkovitz et al. 2001). Since no pre-processing is performed, the statistics are not spoiled and full benefit of statistical modeling is retained. Similarly, since there is no pre-reconstruction arc correction, there is no preprocessing interpolation blurring. Computationally, listmode methods may be more efficient than histogram-based methods for data which are low-count or sparse. On the other hand, histogram-based methods benefit from shared computations of grouped projections and backprojections, and accelerated ordered-subsets implementations are well understood for histogram-based but not listmode algorithms.

In this work we present a histogram-based reconstruction framework that (ideally) achieves exact maximum-likelihood estimation, that shares many characteristics and advantages of listmode reconstruction methods, and that retains the computational advantages of histogram-based reconstruction. In the following sections we describe the proposed LOR-ML reconstruction framework, discuss issues related to comprehensively incorporating all measurement effects into the reconstruction matrix, implement a practical LOR-OSEM algorithm for 2D PET using a rotation-based projector, and evaluate the performance of LOR-OSEM versus conventional OSEM and AW-OSEM in terms of several basic image quality measures.


2.1. LOR-ML Reconstruction Framework

Full utilization of the ML statistical model requires that the data being reconstructed retain exact Poisson statistics, hence all corrections for system measurement effects need to be incorporated into the reconstruction itself rather than being applied as pre-corrections. Figure 2 shows data processing schemes for conventional MLEM, AW-MLEM, and the proposed LOR-MLEM algorithm (and corresponding OSEM versions). The LOR-based approach preserves exact Poisson statistics and reduces interpolation error by combining the arc correction and projection operations. The LOR-ML algorithm is very closely related to listmode ML, and in some cases they are mathematically identical as shown in the Appendix. The difference between LOR-MLEM and conventional MLEM lies in a subtle but important difference in the definitions of the measured data [p with tilde] and the system matrix F.

In order to reconstruct directly from raw LOR histograms, all aspects of the measurement process need to be modeled/incorporated into the reconstruction matrix. This includes effects such as deadtime, detector sensitivity/non-uniformity, randoms, attenuation, and scatter. Furthermore, LOR-ML specifically requires projection to the actual geometry of the detector (as opposed to an evenly-spaced projection process). In this work we do not propose new models for any of these processes; rather, we implement existing compensation methods into our reconstruction matrix rather than pre-correct for these effects.

Along the lines of Vandenberghe et al. (2002), measurements effects are separated into three classes: multiplicative effects, additive effects, and geometric effects. Multiplicative effects include deadtime, sensitivity and non-uniform attenuation, and result in a simple scaling or sensitivity change for individual LORs. The form of the update equation allows these effects to be factored out of the projection/backprojection operations, and they only need remain explicitly in the normalization term. Additive effects such as scatter and randoms introduce broad, low frequency backgrounds which contain very poorly conditioned tomographic information. They can be compensated for by adding a priori estimates to the forward projection (Lange and Carson 1984, Mumcuoglu et al. 1996), leading to very-nearly (but not quite) the ML solution. Finally, geometric effects include the projection formation process or Radon transform as measured by the particular geometry of the detector, plus resolution related effects such as positron range, non-collinearity, and depth-of-interaction (Cherry et al. 2003). Geometric effects should be modeled explicitly in F, and approximations or inaccuracies in such models may give rise to degraded resolution in the reconstructed image. Equation 2 provides the update equation for LOR-OSEM with explicit representation of how the multiplicative, additive, and geometric effects are incorporated:


where Gji represents the geometric components of F, kl is the combined multiplicative factor for LOR (projection histogram bin) l, and âj is the estimated contribution of additive effects to projection LOR j. Note that âj should be the expected value (not noisy estimate) of the summed additive effects to LOR j. If noisy estimates are used, such as estimating randoms using a delayed coincidence window without smoothing, then the additive effect should be handled in a somewhat different manner.

2.2. Implementation Using a Rotation-Based Projector

Using the above formulation, G can either be pre-calculated or computed on-the-fly using a projector (and backprojector). We’ve implemented a rotation-based projector for ring PET tomographs which maps from a regular voxelized image matrix to unevenly-spaced LORs. The projector rotates the image matrix about the z-axis to angle [var phi] using the fast 3-pass method of shears (Di Bella et al. 1996, Paeth 1986), with resampling to unevenly-spaced LORs during the third shear as shown in figure 3. Projection is completed by summing the columns of the rotated matrix (which are now exactly aligned with the raw LORs). The ability to incorporate the resampling directly within the third shear of the rotator offers computational advantages and reduces interpolation errors as compared to performing a separate arc correction. Other projectors, such as ray- or line-driven, could likewise map directly to unevenly-spaced LORs. We chose the rotation-based projector for its computational efficiency, and also because it provides a convenient framework for image-based modeling of depth-dependent geometric effects (e.g. non-collinearity) and scatter.

Figure 3
The rotation-based projector can combine the resampling to unevenly-spaced LORs with the third pass of the 3-pass method-of-shears rotator. In the conventional rotator, the third pass involves a row-by-row shift in x. The LOR-based version combines the ...

Advanced Issues

When imaging on a ring tomograph with an even number of detectors, the number of LORs for projection angles [var phi] alternate being odd and even depending on whether or not the centermost crystals are exactly centered (i.e. at s=0), or are offset by half a bin (Defrise and Kinahan 1998). The odd and even bins from successive angles are often combined in a process called interleaving, generally with the approximation that they all come from the same angle. When using the rotation-based LOR projector, no approximation is necessary as the exact positions of each LOR can be used during the resample-to-LOR-spacing step. Likewise, when LORs with δ ≠ 0 are included in 2D data (i.e. span > 1), the resultant depth-dependent axial blurring can be modeled in the projector using a 1D convolution kernel in a manner analogous to collimator-response modeling in SPECT.


The LOR-OSEM algorithm was evaluated using a series of simulations and phantom experiments designed to characterize the spatial resolution, contrast, and noise properties of the reconstructed images. The LOR-OSEM results were compared with both conventional pre-corrected OSEM and conventional AW-OSEM.

3.1. NCAT Simulation Experiments

The NCAT torso phantom (Segars et al. 1999) was used to create a digital representation of an 18F-fluorodeoxyglucose (FDG) activity distribution in a medium-sized torso typical for whole-body oncologic imaging. Activity concentration in the liver: myocardium: lung: bone: soft tissue was set to 2.8:1.3:0.4:1.6:1.0. A 35.2 cm (transverse axis) × 23.3 cm (anterior-posterior axis) × 41.7 cm (torso height) phantom was created using cubic 0.0493 cm “sub-voxels”, and then collapsed to form 256 × 256 transaxial slices with a pixel size of 0.197 cm and a slice thickness of 0.394 cm. The rotation-based projector was used to simulate 2D PET LOR measurements of the chest for a ring tomograph matching the geometry of the Advance scanner (General Electric Medical Systems), producing a 283 LOR × 35 slice × 336 angle histogram. The effects of non-uniform attenuation were modeled, and scatter and randoms were estimated using a simplified convolution-based procedure based upon (Bailey and Meikle 1994). Use of the same projector for data simulation and reconstruction isolated differences between pre-reconstruction radial-positioning and LOR-based reconstruction, allowing us to directly study this effect. Differences in spatial resolution in the presence of all geometric factors such as positron range, non-collinearity, and depth-of-interaction were studied experimentally using the hot rod resolution phantom described in the next section.

The LOR histogram projection data were scaled to a total of 6 × 106 total prompts over the 35-slice field-of-view for a single bed position, including 4.95 × 106 trues, 0.75 × 106 scatters, and 0.30 × 106 random coincidences. An ensemble of 20 independent realizations of Poisson noise were simulated using a pseudo-random number generator. Six sets of projection data were prepared for each noise realization and for the noisefree case: (1) raw LOR histogram data with no pre-processing; (2) LOR histogram data pre-corrected for scatter and randoms by subtracting the noise-free scatter and randoms components; (3 and 4) the datasets from (1) and (2) were interpolated to 256 evenly-spaced 1.97 mm projection bins (i.e. arc correction); and (5 and 6) the 256-bin projection matrices were collapsed to provide 128-bin datasets with 3.94 mm bins.

The data were reconstructed using three algorithms: (1) conventional OSEM with pre-correction for attenuation, scatter, and randoms; (2) AW-OSEM, where attenuation, scatter and randoms were incorporated into the reconstruction as in eq. 2; and (3) LOR-OSEM, reconstructing directly from unevenly-spaced LOR histograms with all effects incorporated as in eq. 2. For each algorithm, both 128 × 128 and 256 × 256 images were reconstructed using 3.94 and 1.97 mm voxels, respectively. The images were reconstructed out to 4 iterations with 28 subsets, storing the result for each subiteration. No post-reconstruction filter was applied.


Variance images over the 20 noise realizations were calculated, and noise was characterized as the average s.d./mean for voxels in an 8 cm region surrounding the heart. Spatial resolution was characterized by placing a point source 10 cm off-axis and repeating the noisefree data simulation and reconstructions. The images with and without point source were subtracted to produce point source-only images, and resolution was computed as the full width at half maximum (FWHM) of gaussian fits to radial profiles of the point source. The resolution and noise measures were plotted as a function of iteration/subiteration, demonstrating differences in rates of iterative recovery of spatial resolution and noise. In order to separate convergence rates from the comparison, the results were re-plotted as average s.d./mean vs. spatial resolution, which in some sense may be interpreted as an analysis of the noise versus signal in the image.

3.2. Hot Rod Resolution Phantom

Spatial resolution was evaluated experimentally using the Deluxe Jaszczak Phantom with hotrod insert (Data Spectrum Corp.). The phantom contained six wedges of hot rods with diameters 4.8, 6.4, 7.9, 9.5, 11.1, and 12.7 mm, where the center-to-center rod spacing in each wedge was twice the rod diameter. The phantom was filled with 74 MBq 18F-FDG in water and positioned on the imaging table of the Advance scanner with the rods aligned with the long axis of the scanner. The phantom was centered in the field-of-view and scanned for one hour in 2D mode using twelve 5 min. timeframes. Delayed coincidences were stored in separate files to be used later for randoms compensation. A transmission scan was acquired using rotating 68Ge sources for 10 min. The scan data, plus all scanner calibrations, normalizations, randoms and scatter estimate were then offloaded to a Linux workstation for subsequent offline processing.

The data from the 12 timeframes were added to produce a high count dataset with low noise. The transmission data, delayed coincidences, scanner calibrations and normalizations, and scatter estimate data were processed into multiplicative, additive, and geometric components as described in section II.B. The data were then reconstructed onto both 128 × 128 and 256 × 256 image matrices using 2 iterations 28 subsets of conventional OSEM, AW-OSEM, and LOR-OSEM in the same manner as described for the simulation experiments. No post-reconstruction filter was applied. The ten centermost slices of the hot rod portion of the phantom were summed to provide low noise images, and profiles across the smallest rods were drawn to facilitate comparison.

3.3. Whole-Body Tumor Phantom

The AW-OSEM and LOR-OSEM algorithms were compared for reconstructing whole-body tumor images using a whole-body phantom with 22Na lesions (Kadrmas and Christian 2002). The phantom consisted of a 3D Hoffman brain phantom (Data Spectrum Corp.); an anthropomorphic thorax phantom with liver, lungs, and realistic skeletal attenuation (Radiology Support Devices, Inc.); and an elliptical pelvis phantom with bladder. Twenty-seven 22Na lesions were mounted throughout the thorax of the phantom on monofilament line. The lesions consisted of Lucite spheres machined to inner diameters of 7, 8, 12, and 16 mm (measured volumes of 0.17, 0.27, 0.91, and 2.10 mL) filled with 22Na epoxy mixtures. Four activity concentrations were used for each size: 5.55, 8.36, 13.91, and 22.24 kBq/mL on the day of the experiment. The phantom background compartments were filled with 81 MBq of 18F-FDG solutions in water to produce concentrations relative to soft tissue of 6.6:1 in the brain (gray: white = 4:1), 2.0:1 in the liver, 0.37:1 in the lungs, and 23:1 in the bladder. This setup provided an activity distribution representative of whole-body FDG oncologic imaging, including lesions of four different sizes and numerous target: background ratios ranging from 1.8:1 to 25:1 in the various phantom compartments.

The phantom was imaged in 2D mode on the Advance scanner using six bed positions, each including a 7 min. emission and 3 min. transmission acquisition. Delayed coincidences were stored in a separate file, and all data including scanner calibrations and normalizations were offloaded for subsequent processing. The data were reconstructed using 2 iterations 28 subsets AW-OSEM and LOR-OSEM as described earlier, and the reconstructed images were smoothed with a 4.6 mm Gaussian post-filter. Reconstructed background levels in the images for the two algorithms were compared and found to be identical; given this, lesion contrasts were compared for the two algorithms by computing the peak voxel value for each lesion in the reconstructed images. Gaussians were also fit to horizontal and vertical profiles of each lesion, and the fitted FWHM for each lesion was computed as a figure-of-merit closely related to spatial resolution and partial-volume effect.


4.1. Simulation Results

Figure 4 shows example images from the phantom simulation study, including noisefree and noisy images, an ensemble variance image, and a s.d./mean image. Quantitative analysis of these results are shown in figure 5, which plots the spatial resolution and noise measures as functions of iteration (subiteration) and of each other. The spatial resolution results indicate that LOR-OSEM recovered better spatial resolution than AW-OSEM and conventional OSEM, and also that spatial resolution is recovered with iteration at different rates for each of the algorithms. This confirms our hypothesis that LOR-OSEM offers better spatial resolution by reducing interpolation errors associated with arc correction. Recall that the simulation study was setup to isolate this component of spatial resolution, and that components of spatial resolution unrelated to the differences between LOR-OSEM and the conventional OSEM algorithms were not included in the simulations.

Figure 4
Example images from the simulation phantom experiment: (top row, left to right) noisefree LOR-OSEM reconstruction with point source, without point source, and point source-only difference image; (bottom row, left to right) example single noise realization, ...
Figure 5
Results of the phantom simulation analysis of spatial resolution and noise for 128 × 128 images (left column) and 256 × 256 images (right column). The top two rows show plots of the spatial resolution and noise measures, respectively, ...

The plots on the second row of figure 5 show the noise measure plotted as a function of iteration for each algorithm. Care should be taken in evaluating these plots since each algorithm had a somewhat different rate of iterative convergence. A better comparison can be made by studying the bottom row of figure 5, which shows noise versus spatial resolution. These plots show that LOR-OSEM provided roughly the same degree of improvement over AW-OSEM that AW-OSEM provided over conventional pre-corrected OSEM. This result was due in part to the improved spatial resolution properties of LOR-OSEM, and it is also consistent with the hypothesis that LOR-OSEM provides lower noise levels by retaining exact Poisson statistics in the data and taking full advantage of the ML statistical model.

4.2. Experimental Results

Figure 6 shows images and profiles of the hot rod resolution phantom reconstructed by conventional OSEM, AW-OSEM, and LOR-OSEM. LOR-OSEM provided better resolution of the smaller rods, as indicated by the higher peaks and lower troughs in the profiles. This difference was more pronounced in the 128 × 128 images, where LOR-OSEM by nature reconstructed from the full 283 LORs but the conventional algorithms reconstructed from 128 bins. However, the LOR-OSEM results were also measurably better for the 256 × 256 images where all algorithms reconstructed from a similar number of projection measurements. These experimental data are consistent with the simulation results. Since these results demonstrate consistent differences between conventional OSEM, AW-OSEM, and LOR-OSEM, we focus our remaining analysis upon differences between the currently-accepted AW-OSEM and proposed LOR-OSEM algorithms.

Figure 6
Reconstructed images of the hot-rod resolution phantom and profiles across the smallest (4.8mm) rods for 128 × 128 (top) and 256 × 256 (bottom) images. The profiles show sharper resolution for LOR-OSEM as compared to the conventional algorithms. ...

Example images from the whole-body tumor phantom experiment are shown in figure 7. The images reconstructed with AW-OSEM and LOR-OSEM showed similar overall activity distributions with no significant artifacts, and there were some differences in noise texture and lesion clarity. Differences in reconstructed lesion sizes and contrasts are presented in more detail in figure 8, which compares the peak voxel value and profile FWHM size for each lesion for AW-OSEM versus LOR-OSEM. The results show a trend toward somewhat higher lesion contrast and smaller reconstructed size for LOR-OSEM as compared to AW-OSEM, which again is consistent with the result that LOR-OSEM offers better spatial resolution (hence less partial-volume effect upon reconstructed lesion size).

Figure 7
Example maximum-intensity projections (top) and transaxial slices (bottom) of the whole-body tumor phantom reconstructed with AW-OSEM (left MIP, top row transaxials) and LOR-OSEM (right MIP, bottom row transaxials).
Figure 8
Scatter plots of peak lesion intensity (left) and FWHM size (right) for LOR-OSEM versus AW-OSEM. Since the reconstructed backgrounds for the two algorithms were the same, the lesion height plot correlates with tumor contrast in the images. The LOR-OSEM ...

4.3. Projection and Reconstruction Times

Table I shows projection and reconstruction times for AW-OSEM and LOR-OSEM measured on a 1.8 GHz 64bit AMD Opteron Linux workstation (cost ~2,000 USD). All reconstructions used the full 336 angles and 35 slices inherent to the Advance scanner used throughout this paper. Projection and reconstruction times for LOR-OSEM were longer than for AW-OSEM, primarily because of the increased computation required for projection to and backprojection from the unevenly-spaced LORs. Note that the times provided do not include the pre-processing time for AW-OSEM, hence the difference in total processing times would be somewhat less than those shown in the table. The reconstruction times for both algorithms are less than the time required to scan each bed position, and the time for LOR-OSEM to reconstruct even a 256 × 256 image is well within clinically acceptable limits.

Table I
Projection and Reconstruction CPU Times


The proposed LOR-based reconstruction framework offers many of the same advantages that listmode ML reconstruction offers, and it can be easily applied to existing tomographs that do not offer listmode data access. For discrete detectors where there is a one-to-one correspondence between possible event attributes and LOR histogram bins, the LOR-ML and listmode ML algorithms are mathematically identical. The relative advantages of LOR-ML versus listmode ML lie in differences in data storage efficiency and differences in implementation-specific computational issues. For example, the LOR-ML approach offers potential computational benefits in sharing computations for all LORs at a given angle, and it also carries the advantage of retaining the familiar histogram-based architecture for which subsetting rules and regularizations procedures are well understood. Listmode ML, on the other hand, may offer benefits for storing low count or sparse datasets, and it may be preferable for non-discrete detectors (e.g. rotating tomographs) where binning may result in a loss of data quality.

In comparing LOR-OSEM to AW-OSEM and conventional pre-corrected OSEM, LOR-OSEM was found to provide a better noise-resolution tradeoff than the conventional algorithms. The degree of improvement for LOR-OSEM over AW-OSEM was about the same as for AW-OSEM over pre-corrected OSEM. We attribute this to the combination of the exact statistical model offered by LOR-OSEM and reduction of interpolation errors accompanying arc correction—when arc correction is performed pre-reconstruction, separate interpolation errors cumulate for the arc correction and projection steps; this is reduced to a single set of interpolation errors when the arc correction is combined with the projector. While this work was developed in the context of 2D PET with discrete detectors, many of the principles are applicable to other tomographs, fully-3D PET, and SPECT reconstruction. We conclude that fully modeling the tomograph measurement process during iterative reconstruction, as opposed to pre-correcting for some or all of the measurement degrading effects, offers a more robust statistical reconstruction with favorable noise and resolution properties. Requiring only a change in how existing models for image measurement effects are implemented, LOR-OSEM offers a viable alternative to AW-OSEM which can provide higher quality reconstructions.


The author would like to thank Paul Christian, CNMT for his assistance in acquiring experimental data, Rolf Clackdoyle, Ph.D. for helpful discussions on listmode reconstruction, and Charles Stearns, Ph.D. and Scott Wollenweber, Ph.D. for their assistance in accessing the raw data, normalizations, and calibrations for the Advance scanner. This work was supported by a Benning Seed Grant from the Department of Radiology, University of Utah.

APPENDIX: Proof that LOR-MLEM equals Listmode MLEM for discrete detectors

For a detector with a discrete measurement space, e.g. a non-rotating pixelated PET tomograph operated in static imaging mode with preset upper and lower energy thresholds (so that neither time nor energy are stored for each event), the listmode MLEM update equation can be written:


Here we have taken advantage of the discrete nature of our detector so that there is a one-to-one correspondence between each possible unique set of listmode event attributes and the measurement bins. The index jk identifies the LOR for listmode event k. If we reorganize the backprojection summation by grouping listmode events into sets Ωj = {k|jk = j} for every j=1…M, then we can write:


which is the histogram-based LOR-MLEM update equation. Note that this equality only holds for LOR-MLEM, and that the equality does not hold for conventional algorithms which pre-process the data in any way which alters their Poisson statistical nature. In other words, it is implicit that the same F be (correctly) used for both listmode and histogram-based MLEM in this derivation. Likewise, the equality only holds for a discrete detector where a 1-to-1 correspondence is formed between the listmode event attributes and the histogram bins.


  • Bailey DL, Meikle SR. A convolution-subtraction scatter correction method for 3D PET. Phys Med Biol. 1994;39:411–24. [PubMed]
  • Barrett HH, White T, Parra LC. List-mode likelihood. J Opt Soc Am A. 1997;14:2914–23. [PMC free article] [PubMed]
  • Bouwens L, Van de Walle R, Gifford H, King M, Lemahieu I, Dierckx RA. LMIRA: list-mode iterative reconstruction algorithm for SPECT. IEEE Trans Nucl Sci. 2001;48:1364–70.
  • Browne J, De Pierro AR. A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography. IEEE Trans Med Imag. 1996;15:687–99. [PubMed]
  • Byrne C. Likelihood maximization for list-mode emission tomographic image reconstruction. IEEE Trans Med Imaging. 2001;20:1084–92. [PubMed]
  • Byrne CL. Block-iterative methods for image reconstruction from projections. IEEE Trans Imag Proc. 1996;5:792–4. [PubMed]
  • Cherry SR, Sorenson JA, Phelps ME. Physics in Nuclear Medicine. 3. Philadelphia, PA: Saunders; 2003.
  • Comtat C, Kinahan PE, Defrise M, Townsend DW. Fast reconstruction of 3D PET data with accurate statistical modeling. IEEE Trans Nucl Sci. 1998;45:1083–9.
  • Defrise M, Kinahan PE. The Theory and Practice of 3D PET. the Netherlands: Kluwer Academic Publishers; 1998. Data Acquisition and Image Reconstruction for 3D PET.
  • Di Bella EVR, Barclay AB, Eisner RL, Schafer RW. A comparison of rotation-based methods for iterative reconstruction algorithms. IEEE Trans Nucl Sci. 1996;43:3370–6.
  • Hebert T, Leahy RM. Fast methods for including attenuation correction in the EM algorithm. IEEE Trans Nucl Sci. 1990;37:754–8.
  • Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imag. 1994;13:601–9. [PubMed]
  • Huesman RH, Klein GJ, Moses WW, Qi J, Reutter BW, Virador PR. List-mode maximum-likelihood reconstruction applied to positron emission mammography (PEM) with irregular sampling. IEEE Trans Med Imag. 2000;19:532–7. [PubMed]
  • Kadrmas DJ. Fully-3D PET reconstruction for raw LOR histograms (abstract) J Nucl Med. 2003;44:163. [PubMed]
  • Kadrmas DJ, Christian PE. Comparative evaluation of lesion detectability for 6 PET imaging platforms using a highly reproducible whole-body phantom with (22)Na lesions and localization ROC analysis. J Nucl Med. 2002;43:1545–54. [PubMed]
  • Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomogr. 1984;8:306–16. [PubMed]
  • Levkovitz R, Falikman D, Zibulevsky M, Ben-Tal A, Nemirovski A. The design and implementation of COSEM, an iterative algorithm for fully 3-D listmode data. IEEE Trans Med Imag. 2001;20:633–42. [PubMed]
  • Mumcuoglu EU, Leahy RM, Cherry SR. Bayesian reconstruction of PET images: methodology and performance analysis. Phys Med Biol. 1996;41:1777–807. [PubMed]
  • Paeth A. A fast algorithm for general raster rotation. Graphics Interface. 1986:77–81.
  • Parra L, Barrett HH. List-mode likelihood: EM algorithm and image quality estimation demonstrated on 2-D PET. IEEE Trans Med Imag. 1998;17:228–35. [PMC free article] [PubMed]
  • Reader AJ, Erlandsson K, Flower MA, Ott RJ. Fast accurate iterative reconstruction for low-statistics positron volume imaging. Phys Med Biol. 1998;43:835–46. [PubMed]
  • Segars WP, Lalush DS, Tsui BMW. A realistic spline-based dynamic heart phantom. IEEE Trans Nucl Sci. 1999;46:503–6.
  • Shepp LA, Vardi Y. Maximum likelihood estimation for emission tomography. IEEE Trans Med Imag. 1982;1:113–21. [PubMed]