Analytic time-of-flight positron emission tomography reconstruction: three-dimensional case

In a positron emission tomography (PET) scanner, the time-of-flight (TOF) information gives us rough event position along the line-of-response (LOR). Using the TOF information for PET image reconstruction is able to reduce image noise. The state-of-the-art TOF PET image reconstruction uses iterative algorithms. This study introduces an analytic TOF PET algorithm that focuses on three-dimensional (3D) reconstruction. The proposed algorithm is in the form of backprojection filtering, in which the backprojection is performed first by using a time-resolution profile function, and then a 3D filter is applied to the backprojected image. For the list-mode data, the backprojection is carried out in the event-by-event fashion, and the timing resolution determined weighting function is used along the projection LOR. Computer simulations are carried out to verify the feasibility of the proposed algorithm.


Introduction
One of the advantages of using time-of-flight (TOF) technology is its ability to reduce the image noise [1,2], and analytic algorithms are able to reconstruct TOF positron emission tomography (PET) images [1][2][3][4][5][6][7][8][9][10][11]. The state-of-the-art TOF PET image reconstruction methodology is to use the iterative algorithms such as TOF ordered-subset expectation-maximization algorithms [2]. The filtered backprojection (FBP) algorithm is not a preferred method nowadays, due to the concerns of potential larger noise amplification with the FBP algorithm than the iterative algorithms. These concerns are not well-founded. As we demonstrated before, the FBP algorithm should perform as well as an iterative algorithm when the iteration number is emulated, and the projection noise is modeled in the FBP algorithm [12]. We believe that analytical image reconstruction algorithm can achieve the same noise level as a linear iterative image reconstruction algorithm, e.g., the iterative Landweber algorithm. The same can be said to the backprojection filtering (BPF) algorithm, which is an analytic algorithm that performs backprojection first and then performs filtering [13]. For the list-mode data, it is computationally more efficient to use a BPF algorithm than an FBP algorithm. We recommend use of a BPF algorithm so that it is fast, robust and rebinning error free. In the conventional BPF algorithm, the backprojected image does not have a finite support, and this makes the final filtering step not exact. However, for a TOF backprojector, the backprojected image has a finite support if the backprojection weighting function has a finite support. The TOF BPF algorithm has a potential to have better accuracy if the TOF information is used. We will show in the later part of this paper that the TOF modified 'ramp filter' is 'more local' than the conventional ramp filter. Here, 'more local' means that the spatial-domain convolution kernel of the filter rolls-off rate is faster. In the BPF algorithm the ramp filter is often referred to as the ρ-filter; we will use the term 'ramp filter' in this study.
Fully three-dimensional (3D) TOF iterative reconstruction is computationally expensive. When the object is completely measured, rebinning methods are available to convert the 3D measurements into two-dimensional (2D) measurements, so that faster 2D image reconstruction can be performed [14,15]. This study will only focus on direct 3D TOF image reconstruction with an analytic algorithm without using rebinning.
The objective of this study is to develop an analytical image reconstruction algorithm for the TOF PET in 3D mode. This current study develops a "first backproject, and then filter" (BPF) algorithm for the 3D TOF PET. This algorithm is computationally efficient as the conventional FBP algorithm and is able to regulate noise as the iterative reconstruction algorithm.
In conventional tomography, the ramp filter is not local and the backprojected image is not zero outside the image support. As a result, the BPF reconstruction is not as accurate as the FBP reconstruction due to the finite size of the backprojection image array size. On the other hand, for TOF PET, the TOF modified 'ramp filter' is more 'local' than the conventional ramp filter, and the TOF backprojection has a finite support. It is expected that the accuracy of the BPF image is better than the conventional BPF algorithm. In the BPF algorithm the ramp filter is often referred to as the ρ-filter; we will use the term 'ramp filter' in this study.

3D TOF BPF algorithm for the '4π' detectors
Let us first consider a hypothetical spherical PET detector that measures LORs from all possible directions in 3D. In this ideal hypothetical case, the sampling geometry is referred to as '4π', because the surface area of a sphere with radius r equals 4πr 2 . Here, we use the 3D spherical coordinate system. Since the point spread function (psf) and the reconstruction filter are spherical symmetric, it is only necessary to specify their expressions in the radial direction.
Thus, for TOF backprojection, the 3D psf can be modified by a time-resolution induced Gaussian profile function as We now proceed to find the 3D Fourier transform of the psf (2). Since the psf is only a function of the radial direction r, the 3D Fourier transform of the psf can use the spherical coordinates and be evaluated as follows Since the 3D psf is spherically symmetric, its 3D Fourier transform is also spherically symmetric. The integration result of Formula (3) is independent of the values of θ ω and ϕ ω . Without loss of generality, let θ ω = 0 and φ ω = 0. Thus, the 3D Fourier transform of the psf can be reduced from Formula (3) to Using the formula 3.952.6 in ref. [17], The tomographic post filter in the 3D Fourier domain is the reciprocal of Formula (5), and thus When σ is large, the error function approaches to constant 1 and the filter H(ω) in Formula (6) tends to the ramp filter. When σ is small, the error function can be approximated as erfð ffiffi ffi 2 p πωσÞ ≈ 2 ffiffi ffi π p ð ffiffi ffi 2 p πωσÞ .
Therefore, the filter H(ω) in Formula (6) tends to a constant, which implies that no filtering is necessary. If we require that Fig. 2 The definition of the Ω ψ and the arc length γ, which is a part of a great circle Using an approximation of erfðxÞ ≈ e x −e −x e x þe −x , a close approximation of Formula (7) is Some examples of (7) and (8) are shown in Fig. 1 in the first and second columns, respectively.

3D TOF BPF algorithm for ring detectors
For a more practical ring detector as shown in Fig. 2, not every direction θ * of the LOR is measured. Let Ω denote the occupied region by the measured directions θ * on the unit sphere.
We use Ω π/2 to denote the entire unit sphere, which is also known as the '4π' case. For Ω = Ω π/2 , the 3D psf, g(r), is given in Formula (2) and the tomographic recovery filter's transfer function, H(ω), is given by Formula (6). If Ω is not the full sphere but Ω = Ω ψ is a belt as shown in Fig. 2, the psf is no longer spherically symmetric. The 3D spherical coordinate system has three coordinates: (r, θ, ϕ), and its 3D Fourier-domain counterpart using the spherical coordinates: (ω, θ ω , φ ω ).
The psf in this case is circular symmetric (i.e., the psf is not a function of ϕ). The psf for Ω = Ω ψ is only a function of r and θ as well [16]: In the following, we will find the 3D Fourier transform G ψ of the psf g ψ given in Formula (9). Once G ψ is found, the reconstruction filter is obtained as It is known that the spatial domain multiplication corresponds to the Fourier domain convolution. The 3D Fourier transform of g ψ (r, θ) can be evaluated by the 3D convolution in the Fourier domain as and [16] Therefore, we have Formula (14) is universal; it is valid for the '4π' case as well. For the '4π' non-TOF case, γ(θ ω ) = π and For the '4π' TOF case, G π/2 (ω) was given in Formula (15), and we show it again here as Formula (16): For the ring detector TOF case, we are unable to obtain a closed-form for the 3D convolution in Formula (14). We use the relationship between Formulas (16) and (15) to give an approximate closed-form expression for Formula (14). If we replace 'π' in Formula (15) by ' γ(θ ω ) ', the non-TOF case is transformed to the TOF case (12). We now suggest that we replace 'π' in Formula (16) by ' γ(θ ω ) ', and we obtain an approximate closed-form expression for Formula (14) as The general ring-detector TOF reconstruction filter can be expressed as Filters (6) and (18) are the main results of this study. Formula (6) is exact for the '4π' detection geometry, and Formula (18) is approximate for the ring geometry. These two expressions are identical for the '4π' case.

Results
The main result as expressed in Formula (14), which is the 3D tomographic filter applied to the 3D TOF list-mode backprojection in PET. This main result contains an error function erf, which is not an elementary function. Our main result can be closely approximated without using the error function, as shown in Formula (8). Figure 1 shows some comparison results between the exact expression and approximated expression with different values of time resolution σ. It is observed from Fig. 1 that the filter approaches to the ramp filter as the time resolution σ is poor (e.g., σ → ∞). The filter approaches to a constant as the time resolution σ is perfect (e.g., σ → 0), and in this case, the TOF backproject itself can provide the exact reconstruct without the need of a tomographic filter.
A 3D Shepp-Logan phantom [18] was used in computer simulations to verify the feasibility of the proposed algorithm. Four ring-detector sizes were simulated: ψ = π/2, 3π/8, π/4, and π/8, respectively. For each detector geometry, three TOF uncertainty  values were considered: σ = 1, 10, and 100, respectively. As a comparison, non-TOF reconstruction were also carried out for ψ = π/8, π/4, 3π/8, and π/ 2, respectively. Computer simulation results are shown in Figs. 3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18. All images are shown in the same format. The first row shows 3 orthogonal central cuts of the TOF backprojected image: y-z plane, x-z plane, and x-y plane, where the z-axis is the PET gantry axis. The second row shows the 3 orthogonal central cuts of the reconstructed image by the proposed algorithm. The third row shows the 3 orthogonal central cuts of the 3D filter H as expressed by Formula (18). The non-TOF filter is the reciprocal of right-hand-side of Formula (12). The true image is shown in Fig. 19 using the same 3 orthogonal central cuts. All images are 128 × 128 × 128.

Discussion
When the time resolution σ tends to infinity, the TOF backprojection is different from the conventional backprojection. The difference is the normalization. In the TOF backprojection, the value being backprojected is distributed along a straight line according a profile function. The total area underneath this profile function is 1, similar to the situation of a probability density function. On the other hand, in the conventional backprojection, this profile function is a constant 1 and the total area underneath this profile function is infinity. When comparing a TOF analytical algorithm to a non-TOF analytical algorithm, one should pay attention to the normalization factor involved in the backprojectors, which are different for the two types backprojectors. Finally, we discuss how the backprojection profile function and the Gaussian function in the  tomography filter are determined. We can assume that the TOF timing uncertainty can be modeled as a Gaussian probability distribution with a standard deviation of σ 1 . The parameter σ 1 is determined by the PET system we are using. The TOF backprojection profile function can also be assumed to be a Gaussian function with a standard deviation of σ 2 .
The system psf as defined in Formulas (6) and (18) is Gaussian with a standard deviation of σ 3 , which must satisfy σ 3 = σ 1 + σ 2 . In Formulas (6) and (18), the parameter σ is σ 3 . The parameter σ 2 is only used in the TOF backprojector's profile function and can be arbitrarily chosen. It is an interesting special case that σ 2 = 0, in which the TOF backprojector simply backprojects an event to a single point in the image domain. In this interesting special case, the backprojector is much faster than the conventional non-TOF backprojector and we have σ 3 = σ 1 .
In conventional non-TOF tomography, we have

Conclusions
This study derives a 3D TOF BPF formula for the listmode PET data. This formula can lead to an efficient reconstruction algorithm. This new formula is an extension to Colsher's formula [16] that was derived for the non-TOF PET.
The key component in the BPF algorithm is the tomographic filter H. One way to obtain this filter is by numeric evaluation. Since we have a closed-form psf, g, we can numerically evaluate the 3D fast Fourier transform of g to obtain G, and then numerically calculate the tomographic filter H as 1/G. This study aims to find a closed-form expression for the filter H. We are able to obtain an exact pression of H for the  '4π' case, and to obtain an approximate expression for the ring detection geometry. The computer simulations indicate that the approximate filter gives fairly good reconstruction when the ring detection geometry's span-angle ψ is large. The approximation tends to exact when ψ becomes π/2. The noise control for the BPF image reconstruction can be achieved by the post-filtering method as developed in ref. [12]. Our future plans include performing real PET data reconstructions.