Rotate-and-Slant Projector for Fast LOR-Based Fully-3-D Iterative PET Reconstruction

doi:10.1109/TMI.2008.918328

IEEE Trans Med Imaging. Author manuscript; available in PMC 2009 December 28.

Published in final edited form as:

IEEE Trans Med Imaging. 2008 August; 27(8): 1071–1083.

doi: 10.1109/TMI.2008.918328

PMCID: PMC2798574

NIHMSID: NIHMS160125

Rotate-and-Slant Projector for Fast LOR-Based Fully-3-D Iterative PET Reconstruction

Dan J. Kadrmas, Member, IEEE

Dan J. Kadrmas, The author is with the Department of Radiology, Utah Center for Advanced Imaging Research (UCAIR), 729 Arapeen Drive, University of Utah, Salt Lake City, UT 84108 USA (Email: kadrmas/at/ucair.med.utah.edu)

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

Abstract

One of the greatest challenges facing iterative fully-3-D positron emission tomography (PET) reconstruction is the issue of long reconstruction times due to the large number of measurements for 3-D mode as compared to 2-D mode. A rotate-and-slant projector has been developed that takes advantage of symmetries in the geometry to compute volumetric projections to multiple oblique sinograms in a computationally efficient manner. It is based upon the 2-D rotation-based projector using the three-pass method of shears, and it conserves the 2-D rotator computations for multiple projections to each oblique sinogram set. The projector is equally applicable to both conventional evenly-spaced projections and unevenly-spaced line-of-response (LOR) data. The LOR-based version models the location and orientation of the individual LORs (i.e., the arc-correction), providing an ordinary Poisson reconstruction framework. The projector was implemented in C with several optimizations for speed, exploiting the vertical symmetry of the oblique projection process, depth compression, and array indexing schemes which maximize serial memory access. The new projector was evaluated and compared to ray-driven and distance-driven projectors using both analytical and experimental phantoms, and fully-3-D iterative reconstructions with each projector were also compared to Fourier rebinning with 2-D iterative reconstruction. In terms of spatial resolution, contrast, and background noise measures, 3-D LOR-based iterative reconstruction with the rotate-and-slant projector performed as well as or better than the other methods. Total processing times, measured on a single cpu Linux workstation, were ~ 10× faster for the rotate-and-slant projector than for he other 3-D projectors studied. The new projector provided four iterations fully-3-D ordered-subsets reconstruction in as little as 15 s—approximately the same time as FORE + 2-D reconstruction. We conclude that the rotate-and-slant projector is a viable option for fully-3-D PET, offering quality statistical reconstruction in times only marginally slower than 2-D or rebinning methods.

Index Terms: Fully-3-D imaging, image reconstruction, positron emission tomography (PET)

I. Introduction

ITERATIVE statistical reconstruction techniques in positron emission tomography (PET) provide a robust framework for modeling the statistical properties of the measurement, modeling the image acquisition process, and incorporating prior knowledge (when desired) about the reconstruction solution. While iterative reconstruction has largely become the standard for PET, complete utilization of such methods has been limited due to large computational demands. This is especially true for fully-3-D PET, where both direct and oblique coincidence lines are measured to produce high sensitivity—but very large—datasets. A variety of approaches have been proposed for iterative reconstruction of fully-3-D PET data. The most fundamental, and potentially highest quality, implementation requires that the raw data be operated upon directly by the reconstruction algorithm, thereby making full use of Poisson-based statistical models and avoiding any unnecessary degradation or blurring accompanying data preprocessing steps (e.g., arc-correction). However, these implementations tend to be the most computationally demanding as well. At the other end of the spectrum, preprocessing methods—such as rebinning fully-3-D data into a set of 2-D sinograms followed by 2-D iterative reconstruction [1]–[4]—can be used to greatly speed the reconstruction; however, such methods tend to spoil the Poisson statistics of the data and/or introduce undesired blurring or other degradations.

The objective of this work was to develop a fast and accurate projector for iterative fully-3-D PET reconstruction that offers the potential for full utilization of iterative statistical reconstruction algorithms, that has low computational cost, and that is amenable to modeling the physics and spatially variant resolution effects of the PET acquisition. A key requirement for the projector (and backprojector) is that it can map directly from the image to the raw coincidence line-of-response (LOR¹) measurement space (and vice versa), such that the projection operation explicitly models the nonuniform spacing of the coincidence lines acquired by modern PET tomographs [5]. We use the term “LOR-based” reconstruction to set apart methods that directly reconstruct the raw coincidence-LORs of the scanner from those that operate on arc-corrected or otherwise preprocessed sinograms. In effect, the “arc correction” is incorporated into the projector/backprojector of the reconstruction itself, avoiding interpolations and other degradations that occur when preprocessing projection data prior to reconstruction. The broad class of LOR-based reconstruction methods includes both algorithms that reconstruct from raw LOR histograms as well as listmode algorithms that operate on an event-by-event basis, and they have been studied by us [6] and other groups e.g., [7]–[9]. By reconstructing directly from the raw LOR measurements, the Poisson statistical nature of the PET data is preserved, and the full benefits of maximum-likelihood (ML) statistical reconstruction [10], [11] can be realized.

The proposed projector, which we name the rotate-and-slant projector, takes advantage of certain redundancies in the fully-3-D measurement geometry to provide very efficient projection to multiple oblique sinogram sets. It is an extension of the 2-D rotation-based projector for parallel-beam projections [6], [12], [13], which essentially resamples the image via rotation so that the columns align with the projection bins at a given azimuthal angle ϕ. The new projector also implements a new LOR sampling scheme which provides volumetric integration over the geometric LOR volumes, and that is readily adaptable to future incorporation of spatially distributed sensitive volumes (i.e., the spatially variant point spread function). In this paper, we present an overview of the projector, describe how models for various object- and system-dependent responses can be incorporated, discuss a number of computational optimizations, characterize their effect upon projection and reconstruction times, and compare the new projector to several existing ones in terms of reconstruction time and measures characterizing features of image quality. A full mathematical description of the projector is also provided as an Appendix.

II. Projector Description

A. Rotate-and-Slant Projector

Fig. 1 shows the (s, ϕ, z, δ) coordinate system used to parameterize fully-3-D LORs for a generic ring PET tomograph. In order to develop a computationally-efficient fully-3-D projector, similarities of the measurement geometry were analyzed to identify projection computations which could be grouped and shared among LORs, leading to the rotate-and-slant projector. Like the rotation-based projector for 2-D tomography, the rotate-and-slant projector works by resampling the image so that the columns of the image matrix align with the projection lines. Projection is then accomplished by summing the columns of the resampled image matrix. In 2-D, this resampling amounts to rotating the image matrix to azimuthal angle ϕ. We use the three-pass method of shears [12], [14] for this rotation, which breaks down the 2-D rotation process into a series of three computationally efficient 1-D shears.

Fig. 1

Coordinate system (s, ϕ, z, δ) used to parameterize an LOR for a generic ring PET tomograph. The z coordinate describes the axial position of the midpoint of the LOR, which falls at the point of closest approach to the central axis of (more ...)

For fully-3-D PET, the 2-D rotation-based projector can be used for projecting to direct sinograms (ring difference δ = 0). For oblique sinograms (δ ≠ 0), the projection rays are not perpendicular to the axis of symmetry of the scanner and lie at some polar angle θ. Here, resampling the image matrix so that the columns align with oblique projection rays requires an axial slant (shear) operation, as shown in Fig. 2. This slant is a 1-D depth-dependent resampling operation that can be performed much faster than the 2-D rotation operation of Step 1. Furthermore, since the image can be rotated once, stored in memory, then slanted repeatedly to all ring differences, the computational cost of projecting to multiple ring differences is a relatively small increase over projecting to a single ring difference. As a result, the rotate-and-slant projector has high computational efficiency for projecting to datasets with large numbers of ring differences. Backprojection, which is the transpose of these operations, is likewise computationally efficient since 1-D slants from all ring differences can be summed prior to a single rotation back to the object space. Projection to all segments is additionally speeded through the use of depth compression, as described in Section II-C. The most important elements of the projector are described below, and a full mathematical description of the projector is provided in the Appendix. Note that the projector computes nearly-exact geometric volume integrals of each individual LOR, as described below and in the Appendix, where departures from the exact geometry are noted and all approximations are small relative to the size of the (much broader) point response function.

Fig. 2

The proposed rotate-and-slant projector follows three major steps. Note that Step 2, slanting to ring difference δ, involves only a 1-D slant, which is very fast and can be done repeatedly to all ring differences after a single rotation from Step (more ...)

1) Volumetric Projection

The interpolation method used to performed the rotation and slants has implications for the type of projection integrals that are computed. Since each pair of crystals in coincidence for a typical ring tomograph samples a spatial region with finite volume, we consider the associated LOR to be volumetric. The projector, then, should integrate this entire volume, rather than a line through it. From a purely geometric standpoint, the volume contained between a pair of crystals in coincidence is a parallelepiped formed by the front faces of the crystals. The actual sensitive volume would be spatially variant and extend beyond this parallelepiped region due to positron range, acollinearity, and depth-of-interaction effects. We treat the sensitive volume of an LOR as being separable into two components: the (uniform-sensitivity) geometric component computed directly by the projector, and a spatially variant component which takes into account the remaining features and can be considered the point response function. Future incorporation of models for the spatially variant point response function could be implemented in the rotate-and-slant projector using image-space convolutions [15] in a manner similar to that used for depth-dependent collimator-detector response modeling in single photon emission computed tomography (SPECT) (e.g., [16]). In this work, we consider only the geometric component.

The geometric response function for an opposing pair of crystals of width d is triangular with full-width at half-maximum (FWHM) of d/2 [17]. The geometric component of the response becomes narrower closer to the side of the ring in the same way that adjacent LORs get closer together with larger |s|. The LOR width used here was thus set to half the edge-to-edge crystal width in the transverse direction, and half the ring spacing in the axial direction, where the width in the transverse direction decreased with larger |s| in the same manner that the LOR spacing decreases for a ring tomograph. This does not translate into improved spatial resolution for larger |s|, however—depth-of-interaction effects make the spatially variant portion of the response significantly wider at larger |s|, and the spatial resolution of the tomograph hence worsens away from the center of the scanner. Also note that tomographs typically “interleave” LORs from immediately adjacent angles as shown in Fig. 1 (right), effectively binning the data with double the radial sampling rate but half the angular sampling rate. We use the conventional approximation that the interleaved LORs all fall at the same angle, though the projector could easily be adapted for data binned without LOR interleaving as well (in which case, successive angular bins would alternate between having even and odd numbers of LORs).

In order to integrate the volumetric LORs, the resampling steps of the rotator and slants are performed using volume-of-overlap calculations. The volume-of-overlap is determined by computing the length-of-overlap in the transverse plane after rotation (equivalent to the transverse area-of-overlap, since the rotated columns are parallel to the LORs and the height of each row is 1.0), followed by the area-of-overlap in the axial plane, the product of which provides the volume-of-overlap. This is closely related to distance-driven projection [18] (which computes the transverse length of overlap between the LOR and each image row—almost an area-of-overlap computation), extended to include a more complete volumetric computation and accounting for the position and size of each individual LOR. Consider, for example, 2-D rotation using the three-pass method of shears for projection to a set of direct LORs (δ = 0) at azimuthal angle ϕ. The objective is to resample the image matrix via rotation such that the columns of the matrix align with the LORs, where the left and right edges of each rotated column coincide with the left and right edges, respectively, of the volumetric LORs. As seen in Fig. 3, the first two shears map square pixels to square pixels of the same size, hence area-of-overlap and linear interpolation are equivalent (see the Appendix for more details). However, in the third step, each row of the image matrix is translated to accomplish the shear and at the same time resampled to align with the LOR volumes. If one were projecting to arc-corrected, evenly-spaced projection bins the same size as the image pixels, this third step would amount to a linear interpolation. However, for projection to raw volumetric LORs (whose width and spacing vary with s), the area-of-overlap is computed for the resampling. Fig. 3 also depicts the use of a depth compression factor for acceleration, as described in Section II-C.

Fig. 3

Illustration of rotation by the three-pass method of shears, where arc-correction and a depth compression factor of 2.0 were used. Projection to unevenly-spaced LORs is efficiently accomplished by incorporating the LOR-resampling directly into the third (more ...)

B. Incorporation of Object- and System-Dependent Effects

The rotate-and-slant projector provides a convenient image-based framework for modeling object- and system-dependent measurement effects. These effects can be separated into three distinct classes [19], [20]: 1) multiplicative effects such as attenuation, normalization, and deadtime, which affect the sensitivity of individual LORs, 2) additive effects such as randoms and scatter, which contribute to the number of prompt coincidences measured but do not provide high-resolution spatial information about the source activity distribution, and 3) geometric effects which map the precise location and geometry of each LOR measured by the tomograph, ultimately including variant point responses caused by positron range, noncollinearity effects, and depth-of-interaction [17]. Incorporation of each of type of effect into the rotate-and-slant projector closely follows the method used for the 2-D rotation-based projector as described in [6], with a few additional issues related to oblique LOR positions as described below. In this work, scanner geometries and standard estimates for multiplicative and additive effects were obtained from the scanners using manufacturer-provided research tools, and these were then implemented offline with our reconstruction code.

1) LOR Position

For fully-3-D data, the exact position of oblique (δ ≠ 0) LORs is somewhat more complicated than for direct (δ = 0) LORs. For example, as shown in Fig. 1, the polar angle θ for a given ring difference δ changes as a function of s for a cylindrical scanner. Similarly, “direct plane” and “cross plane” slices are often stacked by alternating between even and odd ring differences to form a final slice spacing that is one-half the ring spacing. This causes the polar angles for oblique LORs in adjacent slices to alternate in value. The position and orientation of each LOR is easily modeled in the slant step of the rotate-and-slant projector (making it a spatially variant 1-D operation), thereby matching the geometry of the PET tomograph in question and avoiding many of the approximations commonly used in other PET projectors. In common practice, LORs from multiple ring differences are sometimes binned together to reduce the number of segments. This polar undersampling is sometimes referred to as “mashing” or “axial compression.” The resultant segment contains contributions from multiple ring differences, hence the corresponding sensitive volume is depth-dependent (in a manner analogous to single-slice rebinning). In our current implementation, we consider the geometric component of the response (as modeled by our projector) to include contributions from the central ring difference contained in the segment, and leave the depth-dependent component of the response to be treated as part of the spatially variant point spread function (as described earlier in Section II-A).

C. Computational Efficiency and Optimizations

There are a number of symmetries and other numerical matters that can be exploited to greatly improve the computational efficiency of the projector [21]–[23]. We have implemented the rotate-and-slant projector in C and tested several relatively straightforward optimizations as described below, measuring their effects upon projection and reconstruction times using the 3-D ordered-subsets expectation-maximization (OSEM) algorithm. Measurements of processing times were performed on one core of a single-CPU 2.66-GHz Intel Xeon Linux workstation with no parallelization or multithreading, and all times could be significantly shortened using parallel-implementations on multiple processors as desired.

1) Depth Compression

When using the rotate-and-slant projector, the volume-integral computations amount to summing each column of the rotated image. Essentially, each row of the rotated image corresponds to a different “depth” between the detector crystals (due to the cylindrical symmetry of the scanner, the “zero” depth can be considered to fall halfway between the crystals—i.e., at the plane intersecting the axis of the detector cylinder—and depth-dependent effects will often be vertically symmetric about this zero depth). Projection to oblique segments involves an axial slant, which amounts to a depth-dependent z-shift of each row of the rotated image. Since the oblique LORs lie at small angles from vertical, the difference in axial shift for neighboring rows is small, and relatively little depth information is required to perform the slant to very good approximation. We have incorporated a depth compression factor into the rotator which collapses the image vertically (by summing adjacent rows) during the second x-shear, essentially compressing the depth information of the rotated image. As a result, both memory utilization and the number of operations to perform multiple slants to the oblique segments is reduced. This depth compression bears some relation to multilevel 2-D projection/backprojection methods such as described by Brandt et al. [24]. The effects of using depth compression upon projection time and image quality are presented in the results, Section IV-B. Using a very large depth compression factor affects the data in a manner similar to single-slice rebinning (SSRB), whereas a depth compression factor of one exactly retains all depth information.

2) Array Indexing

The orders in which the reconstructed image matrix and projection data matrix arrays are indexed can have a significant impact upon processing time. Accounting for the manner in which the arrays are accessed within the projector routine, the array indices can be ordered in such a way as to access contiguous blocks of memory when possible. For example, the 3-D reconstructed image matrix would typically be indexed with three indices: column i, row j, transaxial slice k. This array would commonly be ordered such that i varies fastest, followed by j, and k would vary the slowest; which we refer to as i : j : k ordering. In other words, the (i, j, k) element of the image array would be at memory location i + N_i × (j + N_j × k), where N_i, N_j, and N_k are the dimensions of the i, j, and k directions, respectively.

The i : j : k ordering is efficient for applications which treat the image one slice at a time, e.g., a transaxial image display program. For fully-3-D projection and reconstruction, however, the image volume is better treated as a whole. Considering the main steps of the rotate-and-slant projector as shown in Fig. 2, several different ordering schemes may be considered. We compared projection and reconstruction times for three ordering schemes: i : j : k, k : i : j, and k : j : i. The degree of improvement for the different ordering schemes was dependent upon the image dimensions. The k : i : j ordering scheme consistently provided projection and reconstruction times that were 41%–66% as long as for i : j : k ordering, with k : j : i falling in between. This was not surprising, given that the innermost projection loop is a 1-D axial resampling procedure, and the k : i : j indexing allowed contiguous memory blocks to be accessed for this loop.

Two different projection data indexing orders were also studied, LOR:slice:angle and slice:LOR:angle. The slice:LOR:angle indexing order was consistently found to be 12%–17% faster, and hence it was used throughout this work. Overall, optimizing the image and projection data array ordering schemes resulted in a factor of ~ 2× savings in projection and reconstruction times.

3) Vertical Projection Symmetry

The approximately cylindrical geometry of most modern PET scanners provides what we refer to as a vertical symmetry for the projection operation. Referring to Fig. 2, Step 2 of the rotate-and-slant projection requires a depth-dependent axial slant of the rotated image matrix. Defining the zero depth to be at the center of the image matrix (i.e., at the depth of the axis of symmetry of the scanner), the slant for rows at positive depths is the reverse of the slant for rows at negative depths. Since the slant requires a row-by-row 1-D linear interpolation, the interpolation factors can be shared for the positive and negative depths. In effect, this replaces two multiplications with two additions for each voxel in the row, and this occurs within the innermost loop of the slanting routine. This resulted in a time savings of about 18%. Other symmetries or optimizations could also be exploited as in [21]–[23], though we believe that those used here provide the most acceleration without unduly increasing the complexity of the projector.

III. Evaluation Methods

The rotate-and-slant projector was configured for the geometries of three common PET tomographs: the General Electric Medical Systems Advance PET scanner (283 LORs × 336 angles, 35 slices) and Discovery ST PET/CT scanner (249 LORs × 210 angles, 47 slices), and the Siemens Medical Solutions TruePoint Biograph PET/CT scanner (335 LORs × 336 angles, 81 slices). Lookup tables containing the edge positions and orientation of each LOR as functions of s, ϕ, z, and δ were computed based on the geometry of each scanner and provided to the projection routines. An analytical experiment and two physical phantom experiments were performed to evaluate the performance of the new projector as described below. No polar mashing (i.e., no “axial compression”) was used for the experimental data, though different mashing factors (and thus different numbers of segments) were studied for the projection and reconstruction times presented in Section IV-D.

A. Projection Accuracy

The projection accuracy of the new projector was tested versus analytically computed volume-integral projections of a modified high-contrast Shepp–Logan phantom consisting of twelve ellipsoids. The results were compared to those for two conventional projectors: 1) a ray-driven projector based on Siddon et al. [25], which computes a projection line integral using the length-of-intersection of the projection ray with each image voxel, and 2) the distance-driven projector of De Man and Basu [18], which takes into account the width of each LOR and uses a length-of-overlap voxel weighting scheme. These projection models represent a progression from a simple 1-D line-integral projection model to a more volumetric model, and lead toward the fully volumetric rotate-and-slant projection model. All three projectors were programmed into the same C code, using precomplier directives to differentiate the relevant pieces of code, and all used the same set of code optimizations when applicable. Projection accuracy was assessed by computing the root mean squared error (RMSE) between the projected values from each projector and analytically-computed volume-integrals for the same LOR geometry.

B. Reconstruction Performance

Fully-3-D iterative reconstruction performance with the proposed and conventional projectors were evaluated using five reconstruction schemes: 1) Fourier rebinning (FORE, [1], [2]) followed by 2-D attenuation-weighted (AW-)OSEM with the ray-driven projector; 2) conventional 3-D AW-OSEM using the ray-driven projector, where the measured projection data were arc-corrected prior to the reconstruction step; and 3-D LOR-OSEM using 3) the ray-driven projector, 4) the distance-driven projector, and 5) the rotate-and-slant projector. The LOR-OSEM methods reconstructed directly from the raw PET histograms with delayed-coincidences stored in a separate file (i.e., used the ordinary Poisson statistical model applied to nonprecorrected data [26]), where all corrections including the arc-correction were incorporated in the projector/backprojector. In all five reconstruction schemes, the exact same scatter+randoms estimates and attenuation × normalization factors were used, though the order of implementation of each differed for schemes 1), 2), and 3)-5) according to the requirements of each. For example, 1) the scatter+randoms estimate was presubtracted prior to FORE for reconstruction scheme; 2) it was arc-corrected and then added to the forward projection of the iterative reconstruction for scheme; and 3)–5) it was added (without arc-correction) to the forward projections of schemes. The different reconstruction schemes were evaluated and compared using experimentally acquired phantom data as described below.

1) NEMA Phantom Experiment

The NEMA image quality phantom was imaged in 3-D mode on a Discovery ST PET/CT scanner. The phantom was filled with a total of 111 MBq ¹⁸F in water, and contained cold spheres of 37- and 28-mm diameter and hot spheres (8:1 target:background) of 17-, 13-, and 10-mm diameter. The phantom was centered in the scanner field-of-view and imaged for 7.5 min, storing delayed coincidences in a separate file for later use. The raw data and scanner normalizations were then offloaded to a Linux workstation and processed offline into multiplicative, additive, and geometric components, as described in Section II-B. The data were then reconstructed onto 128 × 128 × 49 slice images with 3.125 mm in-plane voxels using each of the five reconstruction schemes and analyzed as described below.

2) Resolution and Contrast Phantom Experiment

The Deluxe Jaszczak Phantom (Data Spectrum Corporation, Hillsborough, NC) was scanned on the Advance scanner to characterize hot-object resolution and cold-sphere contrast. The phantom contained a plastic insert of six wedges of hot rods with diameters 4.8, 6.4, 7.9, 9.5, 11.1, and 12.7 mm, arranged with center-to-center spacing in each wedge equal to twice the rod diameter. Six cold spheres of diameters 9.5, 12.7, 15.9, 19.1, 25.4, and 31.8 mm were also placed in the phantom. The phantom background was filled with 50 MBq ¹⁸F-FDG in water, positioned on the imaging table with the rods aligned with the long axis of the scanner, and centered in the field-of-view. High count data were acquired by scanning the phantom for 30 min in 3-D mode with the septa retracted, again saving delayed coincidences in a separate file, and a 10-min transmission scan was acquired using rotating ⁶⁸Ge pin sources. The scan data were then offloaded, processed, and reconstructed onto 128 × 128 × 49 slice images with 3.125-mm in-plane voxels using each of the reconstruction schemes.

3) Analysis Methods

The five reconstruction schemes were comparatively evaluated using parameters describing image spatial resolution, object contrast, and background noise, taking into consideration differences in the rate of iterative convergence resulting from the different LOR models inherent in different projectors. The effect of using depth compression upon axial resolution was analyzed by fitting Gaussians to axial profiles of the 13-mm-diameter hot sphere from the NEMA phantom experiment, where the Gaussian function was selected as an ad hoc function which fit the axial profiles sufficiently well to characterize changes in the width of the response. The FWHM of the fitted Gaussians were analyzed as a function of depth compression factor, and also compared to those for reconstructions using single-slice rebinning. The effect of depth compression on projection accuracy was also evaluated by studying its effect on the RMSE error for projections of the Shepp-Logan phantom.

Three analysis measures were computed for the Deluxe Jaszczak phantom experiment: 1) the average peak/valley ratio for the wedge of 4.8-mm hot rods, which is a measure closely-related to the recovered spatial resolution, 2) the contrast of the 12.7-mm-diameter cold sphere, computed as (BG-FG)/(BG+FG),² yielding a measure where 0.0 and 1.0 reflect no contrast and perfect contrast, respectively, and 3) the standard deviation of 256 spatially separated background voxels, quoted as a percent of the mean value, where the background voxels were randomly spaced over the image background but away from the phantom objects and outer contour. This final measure provides a measure of the statistical noise in the image. These two measures were first analyzed as a function of iteration to identify differences in the rate of iterative recovery of image features for each of the reconstruction schemes. The noise measure was then analyzed as a function of the resolution and contrast measures, providing comparisons which account for the different rates of iterative recovery.

IV. Results

A. Projection Accuracy

Projection accuracy for each of the three projectors was studied using four discretized versions of the Shepp–Logan phantom with different voxel sizes. The phantom was first discretized onto a 1024 × 1024 × 648 array with 0.4-mm voxels in-plane and a slice thickness of 0.25 mm. This was then collapsed in all directions by factors of 2, 4, and 8 to produce phantoms with successively larger voxel sizes and levels of partial-volume effect due to discretization error. Each phantom was projected with the ray-driven, distance-driven, and rotate-and-slant projectors onto 335 LOR × 81 slice × 336 angle projection arrays using the geometry of the TruePoint Biograph scanner. The RMSEs for a significantly oblique segment (ring difference δ = 20) are shown in Table I, quoted as the percent of the mean nonzero projection value.

TABLE I

Projection Accuracy^a(%RMSE) for Three Projectors

For all matrix sizes, the distance-driven and rotate-and-slant projectors had similar accuracy, and the ray-driven projector had larger %RMSE. This was expected since the ray-driven projector uses a line-integral projection model, whereas the other two projectors are more volumetric. Comparing the different columns of Table I provides an assessment of the contribution of projector error and discretization error, which dominates for the larger voxel sizes (e.g., 128 × 128 matrix) and decreases rapidly with decreasing voxel size. The distance-driven projector had slightly better accuracy for the 128 × 128 image, likely because it had slightly smaller interpolation error than the rotate-and-slant projector. On the other hand, the rotate-and-slant projector had slightly better accuracy for the 1024 × 1024 image, where the discretization error was smaller and the truly volumetric nature of the rotate-and-slant projector produced more accurate results than the area-of-overlap based distance-driven projector. Overall, the accuracy of the distance-driven and rotate-and-slant projectors was deemed to be very comparable.

B. Effect of Depth Compression

The use of a depth compression factor with the rotate-and-slant projector reduces the number of computations required for projection to multiple oblique sinograms. However, since the use of such factors results in some loss of depth information, there may be a consequent loss of projection accuracy for oblique LORs which cross different slices at different depths. Like SSRB, depth-compression creates a depth-dependent approximation for oblique LORs. However, while the error for SSRB increases with increasing depth to the edge of the field-of-view, the error for depth-compression only increases as one moves from the center to the edge of each depth-compressed slab (it does not successively increase across all slabs). For example, when using a depth-compression factor of 8 with a 128 × 128 image, we have 16 depth-compressed slabs. The depth-information at the center of each slab is exactly retained, and at the edge of the slab the error is similar to that of SSRB at 8/2 = 4 voxels away from the center of the image. Since depths map (in a sense) to radii in the reconstructed image, the errors map to concentric rings on the reconstructed image as shown in Fig. 4. The error is very small at the central radius of each ring and maximal (though still small for depth-compression factors of ~ 8 or less) between each ring. Note that a very similar analogous approximation exists for FORE, where small ring differences (less than δ_lim as defined in [1]) are binned using the SSRB approximation. A typical value of δ_lim for FORE would be on the order of 4–5, which roughly corresponds to the same approximation encountered for a depth-compression factor of twice δ_lim, or about 8.

Fig. 4

Example images showing the depth-dependent behavior of the approximations introduced by the use of depth-compression. These images show the absolute value percent error of a simulated disc phantom one slice thick, centered in the field-of-view, and reconstructed (more ...)

Fig. 5 shows the CPU times for fully-3-D projection of a 128 × 128 image to LORs for the three scanner geometries (each having a different number of LORs, slices, angles, and segments as listed earlier) with different depth compression factors. Projection times are markedly reduced for depth compression factors up to about 8, beyond which they remain relatively stable. Fig. 5 also shows two figures-of-merit related to projection accuracy—the reconstructed axial FWHM of the 13-mm sphere for the NEMA phantom experiment, and the %RMSE between analytical and projected values for the Shepp-Logan phantom. Both of these measures were largely unaffected by depth compression factors up to 8, beyond which accuracy was degraded due to the loss of depth information. As a reference for comparison, the FWHM of the 13-mm sphere when reconstructed by SSRB followed by 2-D-OSEM was 20.8 mm—somewhat higher that that for the rotate-and-slant projector with the very large depth compression factor of 64. These data suggest that use of a depth compression factor of 8 (for a 128 × 128 image) with the rotate-and-slant projector offers a significant (almost 3 ×) speedup with negligible loss of accuracy. In general, the optimal depth compression factor will depend upon the image dimensions, number of segments, and range of ring differences included in the projection data, and in general would be smaller if very large ring differences are included. All remaining results for the rotate-and-slant projector presented in this paper used a compression factor of 8.

Fig. 5

Effect of using depth compression factors upon projection time (left) and accuracy measures (right). Here, the FWHM characterizes the axial profile of the 13-mm sphere in the NEMA phantom reconstructed image, and %RMSE provides a measure of accuracy for (more ...)

C. Characteristics of Reconstructed Images

Fig. 6 shows example images of the Deluxe Jaszczak Phantom for each of the five reconstruction schemes. The images are shown at 10 iterations OSEM with 21 subsets (16 angles per subset) for each scheme. The 4.8-mm rods are clearly resolved on this dataset, and the smallest 9.5-mm cold sphere is likewise clearly resolved. Small circles of radioactivity surrounding the support posts for the spheres can also be seen between the wedges of hot rods. The images for the rotate-and-slant projector and the other reconstruction methods show similar image quality, and the most significant differences noted were differences in reconstructed noise texture. Horizontal profiles across one row of the 7.9-mm-diameter and 11.1-mm-diameter rods show somewhat better peak-to-valley definition for the distance-driven and rotate-and-slant projectors as compared to the other cases.

Fig. 6

Example images of the Deluxe Jaszczak Phantom at 10 iterations for the five reconstruction methods studied. The smallest 4.8-mm-diameter rods are resolved for each reconstruction method, as is the smallest 9.5-mm-diameter cold sphere. Small circles of (more ...)

Results of the quantitative analysis of the Deluxe Jaszczak phantom experiment are shown in Fig. 7. Plots of the resolution (average peak/valley ratio of the 6.4-mm rods) and contrast (12.7-mm cold sphere) measures versus iteration reveal differences in the rate of iterative recovery of image features for the different reconstruction schemes. The plots on the bottom row of Fig. 7 effectively normalize for this effect, permitting comparison of image noise at the same resolution or contrast. The data demonstrate a trend toward improved image quality measures when moving from the preprocessing reconstruction schemes to fully-3-D LOR-OSEM. This reflects the improved statistical models of the fully-3-D iterative methods, coupled with reduced degradation when arc correction is included in the projector (LOR-based) and more-volumetric projectors are used. It is not clear why FORE outperformed AW-OSEM3D in Fig. 7, though this may have to do with differences in noise correlations which are not well characterized by the noise measure that was used. Notably, the rotate-and-slant projector performed as well as, or slightly better than, the two other projectors studied for the LOR-OSEM 3-D reconstructions. In particular, the rotate-and-slant and the distance-driven projectors provided comparable results for the analysis presented. The rotate-and-slant projector computes the full volume-of-intersection of each LOR with each image voxel, whereas the distance-driven projector accounts for most, but not all, of the volume-of-intersection. However, the rotate-and-slant projector also has a slight amount of blurring associated with the first two steps of the rotation operation which the distance-driven projector does not have. Overall, the accuracy of the rotate-and-slant and distance-driven projectors were deemed very similar, with both projectors providing comparable reconstructions.

Fig. 7

Quantitative analysis of image quality measures for the five reconstruction schemes studied. The resolution and contrast measures (top row) demonstrate broad similarities for each reconstruction method, with some differences in the rate of iterative recovery (more ...)

D. Projection and Reconstruction Times

The CPU times required for fully-3-D projection and 3-D LOR-OSEM reconstruction on a single-CPU Linux workstation with no parallelization or multithreading are shown in Figs. 8 and and9.9. The TruePoint Biograph scanner geometry was used for the projection time computations, and the advance geometry was used for the reconstruction time computations. Both sets of timings were measured as a function of the number of segments included in the reconstruction, where the typical operating points for 2-D and fully-3-D modes are marked. The cases with fewer segments correspond to using larger polar mashing (or axial compression), whereas the cases with the most segments correspond to no polar mashing. It can be seen from Fig. 8 that the projection time for the rotate-and-slant projector includes a component associated with 2-D projection and increases slowly with the number of oblique segments. In contrast, the projection time for the other two projectors rises rapidly with increasing number of projection rays (LORs), and each oblique segment adds significant time to the computation.

Fig. 8

Projection CPU times for the TruePoint Biograph scanner, plotted as a function of the number of oblique segments, for the three projectors studied. A 128 × 128 × 81 slice image matrix was used, and the projectors mapped to 335 unevenly-spaced (more ...)

Fig. 9

Reconstruction CPU times for the Advance scanner, plotted as a function of the number of oblique segments, for 4 iterations OSEM with 14 subsets and four reconstruction schemes. Each sinogram of the raw projection data had 283 unevenly-spaced LORs and (more ...)

The total reconstruction times plotted in Fig. 9 show similar trends. Notably, fully-3-D iterative reconstruction with the rotate-and-slant projector was accomplished nearly as fast as FORE followed by 2-D iterative reconstruction, whereas fully-3-D reconstruction with the conventional projectors took many times longer. Thus, the computational efficiency of the rotate-and-slant projector enables fully-3-D iterative reconstruction to be performed in timeframes approaching those for rebinning followed by 2-D reconstruction. Given the results of the quantitative analysis comparing image quality features, we conclude that the rotate-and-slant projector offers a viable solution for direct fully-3-D PET reconstruction which takes full advantage of statistical algorithms without associated limitations in processing time.

V. Discussion

The proposed projector addresses one of the most important issues for iterative fully-3-D PET reconstruction—that of long reconstruction times. A number of different approaches to accelerating fully-3-D reconstruction have been proposed. The approach most widely used today on commercial scanners consists of rebinning the fully-3-D data to an equivalent 2-D dataset (e.g., via FORE), followed by 2-D iterative reconstruction. While this results in fast processing times, drawbacks include the need for preprocessing which may have some associated blurring or other degradation, along with compromising the Poisson statistical nature of the raw data (though methods of preserving the statistics while rebinning are under investigation). The more fundamentally-direct application of statistical reconstruction to raw LOR data itself (i.e., LOR-based reconstruction) overcomes the drawbacks of rebinning, but has been associated with very long reconstruction times. The proposed projector largely resolves this concern. Other fast fully-3-D projectors such as Fourier-based [27]–[30] and matrix-driven [31] methods are under investigation by several groups, but no direct comparison of projection times or projection accuracy was available for this paper. The rotate-and-slant projector does bear some similarities to the SSP algorithm recently published by Hong et al. [32]. Both approaches rotate the image matrix and exploit projection symmetries to obtain fast projection (and backprojection), though they use very different implementations. Our work explicitly computes volumetric projection integrals but only references the possibility of using parallelization techniques, whereas Hong et al. compute line integral projections but specifically implement single instruction multiple data (SIMD)-based parallelization.

In addition to its computational efficiency and volumetric nature, the rotate-and-slant projector has a number of other desirable properties. It provides a convenient image-based framework for modeling effects such as the spatially variant point spread function and scatter. For time-of-flight (TOF) PET, the TOF depth is immediately available after the rotation and easily accessed. As compared to matrix-driven projectors, it has modest memory requirements and can be used to project and/or reconstruct any voxel size without the need for any special setup (such as the need to recompute a projection matrix for each voxel size). The proposed projector is easily adaptable to a wide variety of scanner geometries and can map directly to unevenly-spaced LORs, something which is more challenging for Fourier-based projectors. It is also well-suited for block-iterative methods, and it offers so-called “embarrassingly simple” parallelization at the angle and/or subset level. Since the rotate-and-slant projector relies upon a rotator, however, there is a small degree of interpolation error associated with the rotation. When reconstructing voxels ~ 4 mm and smaller, these interpolation errors are quite small relative to the size of the PET point spread function, and the analysis results demonstrated accuracy very similar to the distance-driven projector. Perhaps the biggest limitation of the rotate-and-slant projector is that it is not well suited for event-by-event listmode reconstruction because it relies upon sharing operations among grouped sets of LORs.

VI. Conclusion

A computationally efficient projector for fully-3-D PET reconstruction has been developed which can provide fully-3-D iterative reconstruction in times similar to rebinning and 2-D reconstruction, and markedly faster than conventional 3-D projectors. The rotate-and-slant projector was implemented in C and configured for three modern PET scanners, and several optimizations were performed to speed reconstruction times. Several image quality measures were analyzed for experimental phantom data, and the new projector compared favorably with the other projectors studied. The rotate-and-slant projector provides a viable option for performing accurate fully-3-D iterative PET reconstruction with reconstruction times fast enough for routine clinical use. The time savings offered by the new projector reduces the need for aggressive subsetting or use of accelerated algorithms that potentially diminish image quality, and it may also permit the use of more accurate though time-consuming methods of compensating for image-degrading effects.

Acknowledgments

The author would like to thank C. Stearns, S. Wollenweber, and S. Ross at General Electric Medical Systems, and M. Casey at Siemens Medical Solutions, for assistance in understanding the scanner geometries and data formats for each manufacturer. The author would also like to acknowledge E. Frey for previous discussions leading to the indexing schemes used in this work, and thank the anonymous reviewers whose comments and suggestions significantly improved the paper.

This work was supported by the National Cancer Institute under Grant R01 CA107353.

Appendix Mathematical Description of the Projector

Let A_k,i,j represent the N_k × N_i × N_j image matrix value at slice k, column i, and row j. Projection of A to all segments m at azimuthal angle ϕ by the rotate-and-slant projector requires two steps: 1) rotation of A about the z-axis by angle ϕ while simultaneously resampling the rotated image columns to the LOR sampling at angle ϕ and applying depth-compression, and 2) projection to each segment by applying axial shears and summing the resultant image columns. We implement step 1) using the three-pass method of shears as described below, where the LOR-resampling and depth-compression are performed during the 3rd shear. The following description assumes −45° ≤ ϕ < +45°; for angles outside of this range, the image A is rotated as needed by 90° or 180° by swapping and/or inverting the x- and y-axes (i, j indices) so that the remaining rotation angle lies within [−45°, +45°]. Note that other implementations of the rotate-and-slant projector could be used, such as computing the transformation matrix for step 1) and implementing the step as a matrix-vector multiplication. We have not assessed such alternate implementations in this work.

Step 1): Rotation, LOR-Resampling, and Depth Compression

First X-Shear

Let d_y(j) = j − (N_j − 1)/2 be the perpendicular distance of the center of row j from the central axis. The shear is accomplished by shifting each row j (numbered 0, 1,…, N_j − 1)by Δx(j) = −d_y(j) • tan(ϕ/2), using interpolation when Δx(j) is noninteger. Since the shifted and original pixels have the same size and spacing, interpolation by length-of-overlap is identical to linear interpolation, and we have

(1)

where

Here, A′ represents the image matrix after the first x-shear, and floor() returns the operand rounded down to the nearest integer. Equation (1) is applied for every i, j, and k. For efficient computer implementation, the outermost loop would be over j, followed by i, and the innermost loop would be over k.

Y-Shear

The second shear is in the y-direction. Here, d_x(i) = i − (N_i− 1)/2 is the perpendicular distance from the center of column i to the central axis, and Δy(i) = −d_x(i) • sin(ϕ) is the shift for each column i. Length-of-overlap interpolation is again equivalent to linear interpolation, and we have

(2)

where

Here, A″ represents the image matrix after the y-shear. Equation (2) is efficiently implemented with outermost loop over i, middle loop over j, and innermost loop over k.

Second X-Shear + LOR-Resampling + Depth-Compression

The final shear of the rotator is identical to that of (1), but applied to A″ to yield the final rotated image. One could choose to complete the final shear, and then perform the LOR-resampling on the rotated image; however, this would involve consecutive interpolations in the x-direction—one for the final shear, and one for the LOR-resampling. We perform the LOR-resampling simultaneously with the final shear so that a single interpolation is performed, providing a faster implementation with less interpolation error.

For the shear portion of the calculation, d_y(j) = j − (N_j − 1)/2 and Δx(j) = −d_y(j) • tan(ϕ/2) as above. As shown in Figs. 1 and and10,10, the LOR spacing is spatially variant, individual LORs have different widths, and adjacent LORs may be separated by gaps or even overlapping depending on the scanner geometry and LOR model used. In performing the LOR-resampling, we want to compute the intersection of each individual LOR with the rotated image voxels, weighted by the area-of-overlap within the image slice.³ Since the LORs (within the slice) are parallel to the columns of the rotated image matrix, the vertical length-of-overlap between LORs and image rows is always 1.0; however, the length-of-intersection of each LOR with each (shifted) voxel of each row needs to be computed.

Fig. 10

Diagrams showing the intersection of LORs (shaded) with image columns in the transverse plane (left) and with image slices in an axial plane (right), after rotation of the image matrix to angle ϕ. The area-of-overlap between LOR A and image voxel (more ...)

Let

be the positions of the left and right edges of LOR n at angle ϕ expressed in units of image pixels; these would typically be obtained from a lookup table describing the specific scanner geometry. Also, let equation M6

and

be the left and right edges of the shifted image voxels for column i and row j (where the voxel positions have been shifted by Δx(j) for the final x-shear). The rotated, LOR-resampled, depth-compressed image matrix, B_{k,n,j_comp}, is then computed as shown in (3) at the bottom of the page. The first summation performs the depth-compression by summing all rows (depths) j such that floor(j/γ) = j, where γ is the depth-compression factor (γ = 1.2,4,8,…). The second summation is over all i such that shifted image column i intersects LOR_n, and equation M8

is the length-of-overlap (and, equivalently, area-of-overlap) between shifted voxel i and LOR_n in row j in units of pixels. We implement (3) with outermost loop over j, middle loop incrementing i and n in conjunction, and innermost loop over k.

Step 2): Projection of Rotated Image to All Segments

B_{k,n,j_comp} equals the sum of the intersection of image voxels within slice k and depth-compressed slab j_comp with LOR n. Projection to a “direct” segment (i.e., with ring difference δ = 0; equivalently a 2-D projection) is simply accomplished by summing over depths

(4)

where

represents the projected value at slice k for LOR n. Projection to oblique segments with ring difference δ ≠ 0 involves applying an axial shear followed by summing over depths. The choice of axial shear depends on certain interpretations of the segment geometry. Consider, first, the case where the measured data were binned with minimum span; i.e., each segment contains two ring differences. For example, the “first” oblique segment would contain ring differences δ = 2 and δ = 3, where slices in line with the rings of the scanner (sometimes referred to as “direct-plane”) slices have δ = 2, and slices lying halfway between rings (“cross-plane”) have δ = 3. The exact axial shearing operation would differ for direct-planes and cross-planes, which would add a small degree of complexity to the shear computation. In this work, we make the approximation that all slices of the segment fall at the average ring difference (δ = 2.5 in this example). This approximation has no effect at zero depth (center of the scanner), but adds a small degree of error in the shear operation for image rows away from the center. For applications where the highest axial resolution is important the exact shear should be used; however, we note that the axial mispositioning/blurring due to this approximation is considerably smaller than that introduced by binning the data with other than minimal span. Also note that, as described in Section II-B. LOR Position, the polar angle θ for LORs with the same ring difference varies as a function of s, an effect which we include in our scanner geometry lookup tables and account for in the axial shear operation described below.

The second consideration for the axial shear relates to binning oblique sinograms with polar mashing or “axial compression,” where groups of ring differences are binned together in order to reduce the total number of segments. The effect of this is to add a degree of depth-dependent axial blurring to each segment similar to that of single-slice rebinning but smaller in magnitude (according to the mashing factor). In our implementation, we project the image to LORs at the central ring difference of the segment, and suggest that the width of the segment can be considered as a component of the spatially variant point spread function. However, since the rotate-and-slant projector can project to many segments with low computation cost, we suggest that polar mashing generally be avoided and the complete fully-3-D dataset be utilized for reconstruction. Of course, other issues such as memory and disk usage should also be considered when making this determination.

(3)

Projection to oblique segment m, then, requires knowledge of the polar angle θ(m, n, ϕ) for each LOR n at the given azimuthal angle ϕ. This would typically be stored in a lookup table for the particular scanner geometry, where we note that the LOR-dependence arises from the s -dependence of θ for scanners with cylindrical (or similar) geometries. The depth of the center of each depth-compressed slab j_comp is given by d_z(j_comp) = γ × (j_comp − (N_y/γ − l)/2), giving a depth-dependent axial shift of

A typical image model uses contiguous evenly-spaced slices, where the slice thickness equals one-half the ring spacing. The intersection between each oblique LOR and an image row is trapezoidal as shown in Fig. 10. Since θ is in general not large, this trapezoidal intersection spans two slices in almost all cases, and the trapezoidal area-of-intersection can be computed by simple linear interpolation at the depth of the center of the image row. In certain situations when the oblique LOR is centered or nearly-centered upon an image column at the given row, the intersection traverses three slices. In our current implementation, we count the contribution from the two columns nearest the LOR for these cases and ignore the small contributions from the third column; this introduces a small approximation into the volume-of-overlap calculation for these specific cases.

According to the vertical projection symmetry (Section II-C.), the shift for positive depths is the reverse of that for negative depths. Let equation M13

be the index for the depth-compressed slab at the opposite depth of j_comp. Projection to LOR n of slice k, oblique segment m is then computed by

(5)

where

and

We implement (5) with outermost loop over image rows (slabs), middle loop over LORs n, and innermost loop over slices k. Projection to all segments at angle ϕ is efficiently accomplished by repeating (5) for all segments m.

Footnotes

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

¹We use “LOR” to broadly refer to the geometric linkage between a pair of detector elements in coincidence. In increasing order of complexity, an LOR may be approximated/interpreted as a line between the centers of detector elements, as a 2-D strip between elements within a slice, as a 3-D volumetric “tube” connecting the surface areas of two detector elements, or as the spatially-distributed 3-D sensitive volume of image space that has potential to give rise to prompt coincidence events between the pair of detector elements in question.

²FG and BG represent the average voxel values within regions-of-interest drawn within the object foreground and background, respectively.

³The rotate-and-slant projector computes projection integrals weighted by (almost) the volume-of-overlap between each individual LOR and the image voxels. The area-of-overlap in the transaxial plane is computed by this rotator/resampler, and the remaining (axial) component of the volume is computed during the axial shears of step 2)–projection to each segment.

References

1. Defrise M, Kinahan P, Townsend D, Michel C, Sibomana M, Newport D. Exact and approximate rebinning algorithms for 3-D PET data. IEEE Trans Med Imag. 1997 Feb.16(2):145–158. [PubMed]

2. Defrise M, Liu X. A fast rebinning algorithm for 3-D positron emission tomography using John's equation. Inverse Problems. 1999;15(4):1047–1065.

3. Liu X, Comtat C, Michel C, Kinahan P, Defrise M, Townsend D. Comparison of 3-D reconstruction with 3-D-OSEM and with FORE+OSEM for PET. IEEE Trans Med Imag. 2001 Aug;20(8):804–814. [PubMed]

4. Liu X, Defrise M, Michel C, Sibomana M, Comtat C, Kinahan PE, Townsend DW. Exact rebinning methods for three-dimensional positron tomography. IEEE Trans Med Imag. 1999 Aug;18(8):657–664. [PubMed]

5. Defrise M, Kinahan P. Data acquisition and image reconstruction for 3-D PET. In: Bendriem B, Townsend DW, editors. The Theory and Practice of 3-D PET. Norwell, MA: Kluwer; 1998. pp. 11–53.

6. Kadrmas DJ. LOR-OSEM: Statistical PET reconstruction from raw line-of-response histograms. Phys Med Biol. 2004;49:4731–4744. [PMC free article] [PubMed]

7. Barrett HH, White T, Parra LC. List-mode likelihood. J Opt Soc Amer. 1997;14:2914–2923. [PMC free article] [PubMed]

8. Byrne C. Likelihood maximization for list-mode emission tomographic image reconstruction. IEEE Trans Med Imag. 2001 Oct;20(10):1084–1092. [PubMed]

9. Huesman RH, Klein GJ, Moses WW, Qi J, Reutter BW, Virador PRG. List mode maximum likelihood reconstruction applied to positron emission mammography with irregular sampling. IEEE Trans Med Imag. 2000 May;19(5):532–537. [PubMed]

10. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomogr. 1984;8:306–316. [PubMed]

11. Shepp LA, Vardi Y. Maximum likelihood estimation for emission tomography. IEEE Trans Med Imag. 1982;MI-1(2):113–121. [PubMed]

12. Di Bella EVR, Barclay AB, Eisner RL, Schafer RW. A comparison of rotation-based methods for iterative reconstruction algorithms. IEEE Trans Nucl Sci. 1996 Dec;43(6):3370–3376.

13. Frey EC, Ju ZW, Tsui BMW. A fast projector-backprojector pair modeling the asymmetric, spatially varying scatter response function for scatter compensation in SPECT imaging. IEEE Trans Nucl Sci. 1993 Aug;40(4):1192–1197.

14. Paeth A. A fast algorithm for general raster rotation. Proc Graph Interface. 1986:77–81.

15. Engle JWEJ, Kadrmas DJ. Modeling the spatially-variant point spread function in a fast projector for improved fully-3-D PET reconstruction. J Nucl Med. 2007:417P–417P.

16. Tsui BMW, Hu HB, Gilland DR, Gullberg GT. Implementation of simultaneous attenuation and detector response correction in SPECT. IEEE Trans Nucl Sci. 1988 Feb;35(1):778–783.

17. Cherry SR, Sorenson JA, Phelps ME. Physics in Nuclear Medicine. 3rd. Philadelphia, PA: Saunders; 2003.

18. De Man B, Basu S. Distance-driven projection and backprojection in three dimensions. Phys Med Biol. 2004 Jun;49(11):2463–2475. [PubMed]

19. Mumcuoglu EU, Leahy RM, Cherry SR. Bayesian reconstruction of PET images: Quantitative methodology and performance analysis. Phys Med Biol. 1996;41:1777–1807. [PubMed]

20. Vandenberghe S, Asseler YD, van de Walle R, Lemahieu I. 2002 IEEE Int Symp Biomed Imag. Washington, DC: 2002. Correction for external LOR effects in listmode reconstruction for PET; pp. 537–540.

21. Egger ML, Joseph C, Morel C. Incremental beamwise backprojection using geometrical symmetries for 3-D PET reconstruction in a cylindrical scanner geometry. Phys Med Biol. 1998;43:3009–3024. [PubMed]

22. Joseph PM. An improved algorithm for reprojecting rays through pixel images. IEEE Trans Med Imag. 1982 Nov;MI-1(3):192–196. [PubMed]

23. Hong IK, Chung ST, Kim HK, Kim YB, Son YD, Cho ZH. Proc 2006 IEEE Nucl Sci Symp Med Imag Conf. San Diego, CA: 2006. Fast forward projection and backward projection algorithm using SIMD; pp. 3361–3368.

24. Brandt A, Mann J, Brodski M, Galun M. A fast and accurate multilevel inversion of the radon transform. SIAM J Appl Math. 1999;60(2):437–462.

25. Siddon R. Fast calculation of the exact radiological path for a three dimensional CT array. Med Phys. 1986;12:252–255. [PubMed]

26. Politte DG, Snyder DL. Corrections for accidental coincidences and attenuation in maximum-likelihood image reconstruction for positron-emission tomography. IEEE Trans Med Imag. 1991 Mar;10(1):82–89. [PubMed]

27. Brasse D, Kinahan PE, Clackdoyle R, Defrise M, Comtat C, Townsend DW. Fast fully 3-D image reconstruction in PET using planograms. IEEE Trans Med Imag. 2004 Apr;23(4):413–425. [PubMed]

28. Matej S, Fessler JA, Kazantsev IG. Iterative tomographic image reconstruction using fourier-based forward and back-projectors. IEEE Trans Med Imag. 2004 Apr;23(4):401–412. [PubMed]

29. Matej S, Kazantsev IG. Fourier-based reconstruction for fully 3-D PET: Optimization of interpolation parameters. IEEE Trans Med Imag. 2006 Jul;25(7):845–854. [PubMed]

30. Zhang-O'Connor Y, Fessler JA. Fourier-based forward and back-projectors in iterative fan-beam tomographic image reconstruction. IEEE Trans Med Imag. 2006 May;25(5):582–589. [PubMed]

31. Scheins JJ, Boschen F, Herzog H. Analytical calculation of volumes-of-intersection for iterative, fully 3-D PET reconstruction. IEEE Trans Med Imag. 2006 Oct;25(10):1363–1369. [PubMed]

32. Hong IK, Chung ST, Kim HK, Kim YB, Son YD, Cho ZJ. Ultra fast symmetry and simd-based projection-backprojection (SSP) alogrithm for 3-D PET image reconstruction. IEEE Trans Med Imag. 2007 Jun;26(6):789–803. [PubMed]