|Home | About | Journals | Submit | Contact Us | Français|
For modern Time-Of-Flight PET systems, in which the number of possible lines of response and TOF bins is much larger than the number of acquired events, the most appropriate reconstruction approaches are considered to be list-mode methods. However, their shortcomings are relatively high computational costs for reconstruction and for sensitivity matrix calculation. Efficient treatment of TOF data within the proposed DIRECT approach is obtained by 1) angular (azimuthal and co-polar) grouping of TOF events to a set of views as given by the angular sampling requirements for the TOF resolution, and 2) deposition (weighted-histogramming) of these grouped events, and correction data, into a set of “histo-images”, one histo-image per view. The histo-images have the same geometry (voxel grid, size and orientation) as the reconstructed image. The concept is similar to the approach involving binning of the TOF data into angularly sub-sampled histo-projections -projections expanded in the TOF directions. However, unlike binning into histo-projections, the deposition of TOF events directly into the image voxels eliminates the need for tracing and/or interpolation operations during the reconstruction. Together with the performance of reconstruction operations directly in image space, this leads to a very efficient implementation of TOF reconstruction algorithms. Furthermore, the resolution properties are not compromised either, since events are placed into the image elements of the desired size from the beginning. Concepts and efficiency of the proposed data partitioning scheme are demonstrated in this work by using the DIRECT approach in conjunction with the Row-Action Maximum-Likelihood (RAMLA) algorithm.
THE potential benefit of using Time-Of-Flight (TOF) information for Positron Emission Tomography (PET) imaging was recognized as early as 1969, as noted in reviews of PET -, but TOF did not receive much attention until the properties of the fast scintillators CsF and BaF2 were publicized , , indicating their potential use for TOF-PET. Much work was done on TOF-PET in the early 1980s, but some properties of these fast scintillators were far from ideal, and it was soon found that better overall imaging performance could be obtained by non-TOF scanner designs using other scintillators. For many years, TOF-PET received little attention, but interest was revived once again by the introduction of new scintillators with an attractive combination of properties -, including fast timing characteristics, good stopping power, and high light output. There is now strong interest in TOF-PET because of the significant performance improvements that are possible, as shown in evaluations using simulated data (including , ) and experimental measurements on TOF scanners -.
Full realization of the TOF benefits requires proper reconstruction tools. Many important elements of the theory of image reconstruction from TOF-PET data were published more than 25 years ago in two seminal papers , . With the extra data dimension brought by TOF information, typical acquired TOF data become very sparse in (histo-)projection space. Consequently, the most natural TOF reconstruction approaches are list-mode approaches - that accurately (geometrically and statistically) treat each TOF event separately. However, for the typical number of acquired events in a clinical patient study, conventional list-mode reconstructions with accurate modeling are computationally very expensive.
The computational burdens of list-mode approaches are further exacerbated by the calculation requirements for the sensitivity matrix, which has to be calculated specifically for each reconstructed object based on its own attenuation map. The sensitivity matrix is calculated using all possible lines of response (given by detector crystal pairs) passing through the reconstruction region, including those for which no events were detected. For modern PET systems with a large number of possible lines of response (LORs), exact calculation of the sensitivity matrix using straightforward space-based approaches is not feasible for routine use. There are several strategies with various levels of approximation that allow faster computation of the sensitivity matrix, such as using a properly down-sampled LOR space -. The number of needed LOR samples and the level of approximations depend on a number of factors (e.g., data statistics, object characteristics, detector geometry) and have to be carefully chosen so that the errors propagated from the sensitivity image ,  do not negatively impact image quality. Another option involving reasonable approximations, which is used in this work, is a Fourier-based approach. The advantage of such an approach is great computational efficiency without the need to down-sample the LOR space.
To make list-mode reconstruction time clinically feasible (even with a multi-node cluster of computers) various simplifications and/or restrictions are often needed, such as truncating the TOF and/or resolution kernels, limiting the size of the basis functions and of their underlying grid density, limiting the number of iterations, limiting the number of counts and time frames in dynamic, temporal, and motion-synchronized studies, or using simplified data correction approaches. The goal of the proposed approach is to make TOF reconstruction efficient without compromising image quality via such restrictions. The efficiency of the proposed approach - Direct Image Reconstruction for TOF (DIRECT)  - is obtained by directly working in image space. That means that the data and all corrections are placed directly into (histo-image) voxels in image space, and consequently all reconstruction operations are performed in image space. At the same time the proposed approach has the ability to preserve geometrical accuracy (within the limits given by image and angular sampling) similar to list-mode approaches. Several aspects of the proposed approach were inspired by concepts introduced into TOF-PET that appeared in papers from early 1980s. However, our proposed approach involves a new level of use of these concepts. These include placing events directly into image space rather than into projection space , , merging/grouping events into sets of views and tilts  (and recently -), and using maximum-likelihood (ML) based statistical reconstruction , ,  (and recently , , ). We also build on our previous work at the University of Pennsylvania, including the image space reconstruction approach (ISRA) , 3D view-by-view row-action maximum-likelihood algorithm (RAMLA) , Fourier-based reconstruction approaches -, TOF list-mode reconstruction , and TOF reconstruction using histo-projections from limited views , .
The work presented here describes the principles of the DIRECT approach and gives preliminary results based on simulation studies to demonstrate the performance and properties of the proposed technique, compared with list-mode iterative TOF reconstruction.
The DIRECT approach encompasses the data partitioning scheme and the methods needed to apply a reconstruction algorithm to the given data structure. In this section we outline basic concepts and operations for transforming (pre-processing) acquired data into the DIRECT format; in the next section we discuss the operations of a reconstruction algorithm within the DIRECT framework.
The DIRECT partitioning involves two steps: the acquired events are first sorted (“view-grouped”) into a set of views and then histogrammed (“deposited”) into histo-image elements.
In this step, the acquired TOF events are sorted into a set of “views,” where each view compromises a range of LOR directions having azimuthal angles within a specified interval and co-polar angles within a specified interval. Based on the angular sampling requirements, events acquired on systems with good TOF resolution can be sorted into a relatively small number of views, with only a minimal or no loss of spatial resolution , , . Events from each view have a common TOF kernel (same representative orientation and kernel width) and common data correction arrays. Furthermore, unlike in the non-TOF case, for the TOF data any spatial resolution loss due to the angular binning has been shown to be approximately spatially invariant ; this invariance allows efficient modeling of the effect during reconstruction. The view-grouping leads to a markedly less sparse and reduced data space and a reduction in the number of TOF kernels to be considered.
Traditionally, binned TOF events are histogrammed into “histo-projections” (projections extended in the TOF direction through an image space). Histo-projections are geometrically similar to rotated and resampled images as illustrated in Fig. 1-left. In our approach, the acquired events are deposited prior to reconstruction directly into the voxels of the “histo-images,” with one histo-image per view. The histo-images have the same sampling (voxel size) and geometry as the reconstructed image, as illustrated in Fig. 1-right. The same data structures are used for the deposition of the correction factors, such as attenuation, geometric sensitivity, estimated scatter, and random coincidences. Various deposition strategies are possible (e.g., an event is assigned with different weights to neighboring bins), and the effect of the deposition operation can be efficiently modeled in the forward- and back-projectors during reconstruction. The “histogramming” deposition operation based on nearest-neighbor interpolation is similar to the one in the “most-likely-position” reconstruction technique . Deposition using more accurate interpolation kernels better preserves information about the positions of individual events, and this will be the subject of future studies. In this work, we have employed a tri-linear interpolation kernel.
The approach using histo-images has a number of advantages:
The proposed approach and principles allow efficient implementation of both iterative and analytical algorithms; however, in this work we concentrate on iterative TOF reconstruction. The algorithm flowchart shown in Fig. 2 is applicable to most iterative reconstruction algorithms. In this paper, we illustrate the DIRECT approach used in conjunction with the row-action maximum-likelihood algorithm (RAMLA) , ,  with a view-by-view update. All data structures utilized while processing a particular view v (in sub-iteration n) are 3D arrays of the same geometry holding relevant data for view v: x(n) - image array, ωv - multiplicative corrections, kv - forward-projection kernel (including models of TOF and detector/LOR resolution effects, model of the deposition operation, and possibly including other models of 3D image resolution, regularization, and image basis functions, such as blobs ), - forward projection, - data estimate, dv - deposited data, sv - estimate of expected scatter, rv - estimate of expected randoms, - data discrepancy, and - correction image. Note that the actual memory requirement for DIRECT is much smaller than implied by Fig. 2. For example, only three full 3D arrays are needed for RAMLA with view-by-view updating and four arrays for RAMLA using subsets, rather than the seven arrays shown. We further use the following notation: ·, /, and + denote point-wise multiplication, division, and addition between the elements of two 3D arrays, or between a constant and elements of a 3D array, is the 3D convolution operation, and (relaxation parameter) and 1 are scalars. For RAMLA, the individual steps (from Fig. 2, solid lines) for sub-iteration n processing view v are:
The multiplicative correction array ωv includes attenuation and other multiplicative correction factors, such as sensitivity including gaps and axial acceptance, which are treated similarly as in other (projection-based or list-mode) statistical algorithms, as discussed in more detail in Section IV.A. The back-projection operation is the transpose of forward-projection as dictated by statistical iterative algorithms. It is interesting to note that the 3D convolution operations in both forward-projection and in its transpose are equivalent, since in the DIRECT approach the data and image have the same format. General forms of the forward-projection, data estimate, and back-projection operations are the same for most iterative algorithms. On the other hand, exact forms of the discrepancy and update operations are determined by the particular statistical algorithm. A subset-type of strategy will differ only in the update operation, which will include a normalized sum of correction images for all of the views from the given subset. Examples of single view images resulting from the various operations are included in Fig. 3.
Forward- and back-projection operations, which are the computational bottleneck of any reconstruction algorithm, can be efficiently implemented within the DIRECT approach on a standard computer as well as on specialized computer architectures. Spatially invariant resolution kernels can be implemented very efficiently in the Fourier domain via Fourier-based 3D convolutions, as demonstrated in this paper. Spatially variant kernels, however, must be implemented in the spatial domain. Even in this case, very efficient implementation of these operations is possible within DIRECT framework. For example, substantial acceleration of the spatial convolution operations has been reported when using commodity graphics cards .
In typical emission data the true events (having a Poisson character) are distorted and contaminated by a number of physical factors. To make the best use of the acquired data and of our knowledge of the acquisition system, these factors should be included into the reconstruction model. The contamination factors can be divided, by their character and the way they are treated, into multiplicative and additive terms. The multiplicative factors include: attenuation of the annihilation photons by the object, the probability of a pair of crystals detecting an event once they are hit by the photon pair (detector normalization factors), and the geometrical restriction of direction/LORs for which true events are detected (axial acceptance angle, detector gaps). The additive factors include estimates of the expected scattered and random coincidences. Within the DIRECT approach, the correction factors are placed into the same histo-image format as the true events. This format allows not only a very efficient reconstruction process but also efficient generation of the correction factors themselves. For analytical reconstruction approaches the data would have to be pre-corrected using these factors. In statistical reconstruction methods, an attempt is made to preserve the Poisson character of the data as much as possible by including the correction factors inside the reconstruction model through appropriate multiplicative factors and additive terms.
In statistical reconstruction, the sequence of the physical effects that occur as the coincident events are generated and detected is modeled within the system (probability) matrix. The system matrix can be factorized into a sequence of operations  modeling individual stages of this process as discussed in detail in review paper :
where Ppositron models the positron range, Pgeom contains the geometric probabilities (without attenuation) that photon pairs from individual image locations (points in a voxel) reach the front faces of given crystal pairs (LORs), Patt is a diagonal matrix containing attenuation factors, Ptof models the timing accuracy (TOF resolution kernel), Pdet.blur models the accuracy of reporting the true LOR position (referred to as the detector resolution kernel in this work), and Pdet.sens is a diagonal matrix modeling the probability that an event will be reported once the photon pair reaches the detector surface (crystals or gaps) - a unique multiplicative factor for each detector crystal pair (LOR) modeled by normalization coefficients, but including also the detector axial extent and detector gaps. In general, the positron range has a small effect (compared to the other factors) for whole-body scanners, particularly for 18F-labeled tracers. Therefore, we will omit this factor from further discussion in this paper.
In the following, we consider that the emission counts and attenuation coefficients are approximately uniform within the size of a voxel, as given by the system resolution. Under this condition the attenuation factors (in Patt) can be considered to be constant for the events included in the corresponding elements in Pgeom, and consequently, places of Patt and Pgeom can be permuted in the sequence (formal structure of the matrices changes with the permutation). Pgeom, Ptof, and Pdet.blur can then be modeled together. In the following derivations we further consider that after the afore-mentioned process the events are angularly grouped based on the TOF angular sampling requirements. The forward model of the true portion of the view-grouped events (p) is then given by:
where G represents the view-grouping operation, P is the system matrix (as discussed above) and x is the emission image. Efficient implementation of this model within DIRECT framework is derived below using the following symbols (consistent with Fig. 2):
For scanners with a uniform distribution of crystals with no gaps and for data pre-corrected for detector normalization factors (as is often done in practice -), the process of the generation-acquisition-deposition of the true events can be modeled as G P = G PsysPω, where Psys combines the deposition kernel, Pdet.blur, Ptof and Pgeom, and Pω represents the attenuation factors Patt. Note that the following derivation does not preclude adding additional blurring operations to the right, such as an image basis function model or Ppositron, e.g., G P = G PsysPωPpositron. The forward-projection of true events for a particular view v, and voxel m can be written as (details on derivation of the following sequence of expressions are below):
where . This derivation is based on the fact that if the range Nv of directions i used in each view v is less than that dictated by the angular TOF sampling requirement , the image values accessed through the convolution sums over convolution kernel give approximately the same value for any direction i within a given view. Consequently, the convolution sums and kernels can be approximated by a representative convolution over for given view v (second expression). This convolution sum is now independent of i within each view v and the outside sum over the directions i can therefore be moved inside and applied only to the multiplicative correction coefficients prior to the convolution operation (3rd expression). Note that while the local convolution sums over the emission activity x within the resolution kernel regions do not depend on the particular direction within a view, the sum over the attenuation (and other multiplicative) factors does, since the attenuation factors are not local, but are determined by the attenuation on the entire length of an LOR. Therefore the sum must still be performed over all directions within each view v. However, because of the rearrangement of sums allowed by the locality of the emission convolutions, the factors ωv,m can be pre-calculated for each voxel in each view, forming the correction matrix ωv for each view v. Forward-projection for a given view v can then be calculated by point-wise multiplication of the image values by this correction matrix, followed by a forward-projection operation that consists of convolution of the corrected image with the (combined) resolution kernel as shown in the flowchart of Fig. 2 (solid lines).
Attenuation factors on lines i in the projection space can be very efficiently calculated from the attenuation image using Fourier-based forward-projectors , . After combining attenuation and other correction coefficients in projection space, we can calculate the multiplicative correction matrix ω in histo-image format by calculating the sum very efficiently via Fourier-based back-projections; there is a separate back-projection for each view, each using the proper range of directions of correction factors in projection space.
Current commercial scanners have modular detectors with gaps between adjacent modules and with discrete changes of orientation of the detector module faces. A proper model requires the gap information to be included in the system matrix as part of the matrix Pdet.sens, even if the data are pre-corrected for other normalization factors. When there are gaps between detector modules, we cannot drop the factor Pdet.sens from eq. 1 (as was considered for the no-gap case). Pdet.sens is a multiplicative correction matrix, as is Patt, but between these in eq. 1 are the convolution operations representedby Pdet.blur and Ptof. However, efficient implementation of the multiplicative corrections is possible only if all view dependent corrections are combined together in the system matrix sequence. Pdet.sens changes abruptly from LOR to LOR, due to gaps and crystal dependent normalizations, and cannot be moved inside the convolution sums. On the other hand, the attenuation factors change more slowly and typically do not change dramatically within the range given by the detector resolution kernel. This allows us to approximate the system matrix sequence by moving the attenuation correction from before to after the convolution (forward-projection) operations. Note that in conventional iterative reconstruction, the attenuation step is, in fact, routinely done after the forward-projection, which is much more efficient for binned projection data. In the following, we derive DIRECT modeling operations for the case when we combine attenuation coefficients with Pdet.sens after the forward-projection operations (dashed line in the flowchart in Fig. 2): G P = G PωPsys, where Psys combines Pdet.blur, Ptof, Pgeom and can include also the image based resolution kernels , image basis function model, or regularizing functions, and Pω includes attenuation, information on detector gaps, and detector normalization. The forward-projection of the true events for a particular view v and voxel m can then be written (following the same principles as in the previous case) as:
The choice of one of the two multiplicative correction modeling approaches as derived above (eq. 3 vs. eq. 4) will depend on the particular system’s characteristics (e.g., existence of gaps and detector resolution) as well as on the imaging application (speed of changes and/or resolution of the emission and attenuation images) and will be the subject of future studies.
In conventional non-TOF PET, a true coincidence carries the information that an annihilation event occurred somewhere along the LOR connecting the opposing pair of detector elements. In TOF-PET, the information about the position of the annihilation is more localized, since the approximate position along the LOR is also measured. Similarly, in TOF-PET the possible locations of scatter are more localized, since a pair of photons that have undergone one or more scatters can be registered in a TOF “bin” only if the difference between the two path lengths from the annihilation to the detectors is within the timing uncertainty of the true photon pairs contributing to this bin. The single-scatter case is the most important in both conventional (non-TOF) PET and TOF PET , . Much effort has been directed towards the development of efficient model-based methods to estimate the expected contribution of single-scatter events to the measured data -. More recently, a model-based single-scatter estimation algorithm has been developed that takes into account the spatial as well as timing distribution of scattered events for TOF-PET data , , and it is in routine use for clinical imaging at our institution. This algorithm can be modified to calculate the scattered events (sv) directly in histo-image format to be used within DIRECT reconstruction. Care has to be taken that the scatter estimates are treated in the same way as the true event estimates in the forward-projection model, including gap information and consideration of normalized or un-normalized events. Various challenges exist for scatter estimation in general (independent of DIRECT approach), such as modeling of out-of-FOV scatter, but their solutions can be directly imported into the DIRECT (or any other) approach.
Similar to scatter, techniques for randoms estimation are independent of the reconstruction approach, and the estimates of the expected randoms (rv) can be directly deposited into histo-images as is done for the true events. It should be noted that although randoms and scatter estimates are shown as separate arrays in the reconstruction flowchart (Fig. 2), once they are estimated, they can be combined into a common array for more efficient storage.
Full scanner simulations were performed using a Monte Carlo tool based upon the EGS4 simulation package . We simulated a realistic whole body TOF scanner having a 25-cm axial field of view (FOV), ±15° axial acceptance angle, with continuous distribution of crystals on a cylindrical detector surface with no gaps and a 5.8-mm spatial resolution . Simulations were performed for timing resolutions of 300 ps and 600 ps.
The simulated phantoms were 27-cm and 35-cm diameter cylinders having clinically relevant volumes and attenuation factors representative of average and heavy patients, respectively . Simulated data included attenuation effects but not scatter or random events. Hot spherical lesions of size 10, 13, 17 and 22 mm and contrast 4:1 were simulated in two different slices (see Fig. 4). We generated 6 different data sets of true list-mode events, each having 40M counts. Different numbers of counts were then used from the simulated lists: an average count case (20M and 30M counts corresponding to the heavy and average patient cases, respectively), and a very low count case challenging the reconstruction approaches (6M and 9M counts, respectively). For both the average and low count studies, the two count levels represent approximately the same acquisition times for the two patient sizes.
To test the data modeling in the DIRECT approach (as derived in Section IV.A) with respect to abrupt attenuation changes and angular view-grouping, we also simulated a “short phantom” with abrupt axial truncation of attenuation and emission distributions and with hot lesions being located in close proximity to this truncation (see Fig. 4). For comparison purposes we also simulated a 25-cm “long phantom” filling the whole axial FOV of the scanner that was otherwise equivalent to the short phantom.
We used a block version of RAMLA in the DIRECT and list-mode TOF, and in the list-mode non-TOF reconstructions. In all reported experiments the relaxation parameter was set to 1.0 (although RAMLA updates were implicitly relaxed by attenuation and sensitivity matrix values). In the DIRECT approach, we grouped and deposited events into 40 × 3 views - 40 intervals in azimuthal angle and 3 intervals in co-polar angle, except for the study on view-grouping effects where we varied the number of views (10-40 azimuthal and 1-5 copolar intervals). The 40 × 3 set of views was fixed in the rest of the studies based on the view-grouping experiments, so that the view-grouping effects will have minimal impact on the results. Each view represented one subset/block of RAMLA, giving us 120 (geometrically-ordered) updates for one pass through the data in the 40 × 3 view case. Note that the tilted views are not complete because of the limited axial extent of the detector, and consequently the effective amount of update is smaller than for the non-tilted views (centered on co-polar angle zero), or for subsets based on temporal subdivision of data from all directions. In the reported studies we employed TOF kernels representing 300-ps and 600-ps TOF resolutions of simulated data, and modeled a 5.8-mm (FWHM) spatially invariant detector resolution kernel and a 4-mm (FWHM) data deposition kernel. In the comparison studies with list-mode approaches, we also incorporated a model of the blob basis function into DIRECT. The resolution and basis function filters were all applied in the Fourier domain, providing the same computational time independent of the kernel size. In the list-mode approach (with temporal-ordering), the events were subdivided into 60 temporal subsets, each representing one subset/block of RAMLA and each subset covering all views. This number of subsets provided similar convergence speed (contrast and noise increase per iteration) to DIRECT using RAMLA with 120 geometrical subsets and a comparable blob model. List-mode reconstruction used blob basis functions placed on a body-centered cubic grid (radius = 10 mm, shape parameter α = 8.63, grid unit cell size = 8 mm), based on our experience with these particular parameters in our previous studies. The final image array was 144×144×62 with 4 mm3 voxels.
The behavior of the proposed approach in conjunction with iterative reconstruction was investigated using the contrast vs. noise trade-off for a range of iterations and reconstruction parameters. Contrast recovery coefficients (CRC) were calculated for all sphere sizes and locations as: , where ps is the mean value in a 2D circular region of interest (ROI) axially and transversely centered over the sphere (with the ROI having the same diameter as the feature and positioned using a neighborhood search to maximize each sphere’s mean value over all realizations), pb is the mean value in the 2D annular region centered over each sphere (with the inner radius being 7 mm larger than the feature radius and with annular thicknesses being carefully selected to give a similar annular area for any feature size while not being influenced by neighboring features), and c is the ideal contrast value. The transverse centering was done with with sub-voxel accuracy (within 0.5 mm), and the contribution of each voxel value to the ROI or annular region was weighted by the proportion of its area within the ROI or annulus. The average CRC value for a given sphere size was then calculated using the three (eight for 10-mm lesions in slice B) sets of each sphere in each of the six noise realizations, leading to a total of 18 (48) different realizations.
Depending on the study goals, there are several established ways of assessing noise within reconstruction approaches. One of our goals was assessment of effects of the studied approaches on the spatial noise within the images. From qualitative visual inspection of the images and image profiles, the spatial structure of the noise in the DIRECT and list-mode approaches was seen to be very similar (e.g., see Figs. Figs.66 and and7).7). Consequently, we believe that a pixel-to-pixel noise measure captures the noise effects within our studies. In our evaluations, we used features having a range of sizes, and our noise measure over fixed size 50-mm ROI included image variations over the corresponding range of spatial scales. The range of spatial scales of the noise measure was limited from below by the 4-mm voxel size, which is about half of the smallest (10 mm) lesion size. An upper limit on the spatial scale of the noise measure was imposed by setting to 50 mm the diameter of the ROI in which the pixel-to-pixel noise standard deviation was evaluated. This upper limit corresponds to about twice the size of the largest (22 mm) lesion. The noise ROI was centered in the feature slices, and the calculated standard deviation values were averaged over all realizations and normalized by the background mean value. We also evaluated the standard deviation across the realizations. The general conclusions for the two noise assessment approaches were consistent, and we report only on the pixel-to-pixel noise measure.
Fig. 5 shows contrast vs. noise trade-off curves for the various reconstruction approaches studied, using the same algorithm and blob basis function model. It can be seen from the curves that DIRECT with a comparable model is able to match the contrast-noise performance of list-mode TOF reconstruction, but for a much smaller computation time (see Table I). For illustration, we also show the curve for DIRECT reconstruction that additionally includes modeling of the intrinsic resolution of the simulated data (not included in the list-mode reconstructions). The resolution modeling provides an improved contrast vs. noise trade-off, as expected with statistical reconstruction independent of data partitioning approach (whether DIRECT, binned, or list-mode). The same modeling can be incorporated also into the list-mode approach, but at the cost of a several times longer reconstruction time. We emphasize that in the DIRECT approach the reconstruction time is the same independent of the resolution kernel size, since the resolution modeling is done in the Fourier domain. Examples of images and mean central profiles (over six realizations) for the contrast vs. noise levels marked by gray ellipses in Fig. 5 are shown in Figs. Figs.66 and and7.7. DIRECT and list-mode TOF and non-TOF reconstructions with the same model are shown for matched noise levels (bottom three rows). It can be seen that the visual quality of the images and profiles are consistent with observations based on the contrast vs. noise curves. Additionally we show also the image (top row) and profile for the DIRECT reconstruction with resolution modeling at a high iteration number to illustrate that it does not lead to an unstable behavior and does not introduce background artifacts.
In the DIRECT approach, each view-grouped histo-image is considered to have a common TOF kernel. In the simplest case, this kernel is approximated by a single TOF kernel at a fixed orientation in 3D given by the central direction (azimuthal and co-polar angles) of the view (“simple TOF kernel” in the following). More accurate modeling would involve a TOF kernel averaged over the range of directions represented in the individual views. The simple TOF kernel is a good approximation if the angular sampling requirements are fulfilled . For example, for a TOF resolution of 300 ps and a spatial resolution of 4 mm, the minimum required number of azimuthal angular intervals is about 18, and for 600 ps it is about 35. In the co-polar angle it would be approximately 3 and 5 tilts for 300 ps and 600 ps, respectively, for the central plane of our simulated system with a ±15° axial acceptance angle, and fewer away from the center. If the angular sampling requirements are not met, the simple TOF kernel does not accurately model the deposited data, and the angular grouping introduces blurring. However, this blurring is spatially invariant in the TOF case  and can be efficiently incorporated into the reconstruction model. Fig. 8 shows contrast vs. noise trade-off curves for the 10-mm spheres for DIRECT reconstruction using the simple TOF kernel from 600-ps data deposited into variable number of azimuthal (top) and co-polar intervals (bottom). No view-grouping blur was modeled in these experiments, so as the number of views decreases below the sampling requirements (35 azimuthal samples and 5 co-polar samples for the 600-ps case), the contrast vs. noise trade-off starts to deteriorate. The angular TOF sampling requirements are important also for the proper modeling of corrections, as discussed in Section IV.
It is worthwhile to mention here that for algorithms with a view-by-view update such as RAMLA (or the ordered-subsets expectation-maximization algorithm), in which the number of updates per iteration depends on the number of views, using fewer views means that we need to perform more passes (iterations) through the data to get the same number of image update operations. Consequently, those algorithms can achieve similar contrast vs. noise trade-offs for a similar total reconstruction time “independent” (within reason) of the number of views used, as seen from Fig. 8. This in turn allows the use of a sufficiently large number of views without increasing the total reconstruction time. This is unlike the filtered back-projection or maximum-likelihood expectation-maximization (without subsets) algorithms in which the number of views directly affects the total reconstruction time. On the other hand, we do not want to increase the number of views too much either, since a higher number of views also means fewer counts per histo-image and consequently sparser data.
While each individual voxel in the reconstructed image is affected only by the emission activity from a local area as given by the TOF resolution, the attenuation contributions come from the whole LOR length. In Section IV, we showed that even those “long distance” factors can be properly modeled in DIRECT, i.e. it is possible to do the view-grouping and still have accurate modeling. Our studies confirmed that even when we group together events for which the photon pairs travel through a significantly different attenuation, or non-attenuation, regions (as is the case for the events emitted from the 10-mm spheres in the short phantom), DIRECT properly models this situation and the reconstruction quality is not affected. Axially we grouped data into 1, 3, and 5 co-polar intervals and compared performance of reconstructions of the short phantom with abrupt attenuation truncation close to the spheres with those from the long phantom with no truncation. Fig. 9 shows transverse and sagittal images of the short phantom for the reconstruction using 3 co-polar intervals. Axial profiles through the short and corresponding long phantom images are shown in Fig. 10. Both cases provide very similar performance with minimal effects of the object truncation on the reconstruction. Fig. 11 shows CRC curves, as a function of the iteration number, for short and long phantom reconstructions for 1 and 3 co-polar interval cases. The 5 co-polar interval case (not shown) provided very similar comparison and contrast values to the 3 interval case (similar to Fig. 8).
The number of acquired counts, voxel size, and number of views determine the count densities of individual histo-images. For very low count cases it might be reasonable to adapt the size of the image elements or amount of the regularization (changing the resolution-noise trade-off). However, the purpose of this study was to investigate the effect of counts by themselves without doing any adaptation based on count levels. Fig. 12 illustrates contrast vs. noise trade-off curves for different count levels and object sizes. Fig. 13 shows representative images of the central slices for comparable contrast levels (as marked by the shaded bar in Fig. 12). In the lowest count case each histo-image contained on average only 1 count per 5 voxels inside the object. As seen from the images and graphs, the combination of DIRECT and RAMLA approach proved to be robust without any undesirable artifacts even for very low counts.
Fig. 14 illustrates contrast versus noise trade-off curves for different TOF resolutions and object sizes. Fig. 15 shows representative images of individual cases for comparable contrast levels (as marked by the shaded bar in Fig. 14). It is clear from the graphs and visual image quality that better TOF timing resolution improves reconstruction performance, as has been previously shown for list-mode TOF reconstruction .
The computational bottlenecks of any iterative reconstruction are the forward- and back-projection operations, and DIRECT is no exception. Forward- and back-projections are performed in DIRECT using 3D convolutions. For the Fourier-based implementation of those operations, virtually all of the DIRECT reconstruction time is taken by the 3D Fourier transformations - four 3D fast Fourier transforms (FFTs) per each update (or view) of the algorithm, for a total of 480 3D FFTs per iteration in our example. FFT routines can handle this task on the DIRECT structures extremely efficiently even on a standard computer, as seen from the Table I. The memory demands of DIRECT are minimal; the approximate sizes of data structures involved in the studies presented are as follows: 320 MB disk file of list data with 40M events; 300 MB disk file with deposited histo-images (image size 144×144×62, 40×3 views; float values stored in scaled 2B integers), and same for corrections (scatter and random images can be stored also in a down-sampled form, since they are slowly varying); and 20 MB total computer memory (RAM) needed by DIRECT for four 3D arrays (only one view is processed and needs to be in the memory at any time). Even if the image size increases several times, the RAM and disk requirements will be still comfortably within the possibilities of standard computers.
Table I shows representative reconstruction times per iteration of DIRECT and list-mode TOF reconstructions on a 2-GHz Mac G5 single processor (image size 144×144×62). Reconstruction parameters were used as outlined in Section V.B. In the list-mode reconstruction we used blobs on the body-centered grid, no detector resolution modeling was applied, and the TOF kernel was truncated at ±3σ of the Gaussian timing resolution kernel. Reconstruction times of DIRECT do not depend on the number of counts or the sizes of the resolution and TOF kernels. Reconstruction times given for the particular single processor are only illustrative to give insight about the relative computational demands of the approaches. The computational gain of the DIRECT will further increase for higher number of counts; it will also increase for larger image sizes, as expected with the Fourier-based approaches. DIRECT, as well as list-mode approaches, can be efficiently parallelized for a cluster of computers or other parallel computer architectures. DIRECT operations can be very efficiently implemented also on GPUs, but this is outside the scope of the presented work.
The large computational demand of the space-based calculation of the sensitivity matrix results from the fact that all crystal pairs (LORs) must be considered, including those for which no events were detected. For modern PET systems with a large number of LORs, calculation of the sensitivity matrix itself can be computationally more expensive than one iteration through the acquired emission events during the reconstruction stage, as illustrated in Table I. There are several possibilities to speed up the sensitivity matrix calculation, as discussed in the Introduction. Substantial speed-up of the sensitivity matrix calculation within DIRECT geometry has been obtained using the Fourier-based forward- and back-projection operations.
We have presented the Direct Reconstruction for TOF data (DIRECT) approach, working directly in image space and allowing very efficient implementation of reconstruction and data correction operations, without compromising image quality. Our feasibility studies show the potential for about an order of magnitude speed-up in reconstruction compared to our clinical TOF list-mode reconstruction of TOF-PET data for a typical patient study. Furthermore, the DIRECT approach allows modeling of resolution and TOF kernels/filters without any extra computational cost. Our future work will involve comprehensive evaluations using measured data, including scatter and random events, from realistic emission and attenuation objects acquired on our clinical TOF-PET scanner ,  and our research TOF-PET scanner with LaBr3 with improved timing resolution , .
The authors thank the reviewers for their valuable comments and suggestions, and gratefully acknowledge L.M. Popescu for helpful comments on this work and M. Werner for help with practical issues concerning reconstruction software and data.
This work was supported by the National Institutes of Health under grants EB002131 and CA113941.
Samuel Matej, Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104 USA.
Suleman Surti, Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104 USA ; Email: ude.nnepu.dem.liam@itrus.
Shridhar Jayanthi, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA ; Email: ude.hcimu@ihtnayaj.
Margaret E. Daube-Witherspoon, Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104 USA ; Email: ten.nozirev@webuad.
Robert M. Lewitt, Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104 USA ; Email: ude.nnepu.shpu@ttiweL.treboR.
Joel S. Karp, Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104 USA ; Email: ude.nnepu.dem.liam@prakleoj.