An analytic reconstruction method for PET based on cubic splines

PET imaging is an important nuclear medicine modality that measures in vivo distribution of imaging agents labeled with positron-emitting radionuclides. Image reconstruction is an essential component in tomographic medical imaging. In this study, we present the mathematical formulation and an improved numerical implementation of an analytic, 2D, reconstruction method called SRT, Spline Reconstruction Technique. This technique is based on the numerical evaluation of the Hilbert transform of the sinogram via an approximation in terms of 'custom made' cubic splines. It also imposes sinogram thresholding which restricts reconstruction only within object pixels. Furthermore, by utilizing certain symmetries it achieves a reconstruction time similar to that of FBP. We have implemented SRT in the software library called STIR and have evaluated this method using simulated PET data. We present reconstructed images from several phantoms. Sinograms have been generated at various Poison noise levels and 20 realizations of noise have been created at each level. In addition to visual comparisons of the reconstructed images, the contrast has been determined as a function of noise level. Further analysis includes the creation of line profiles when necessary, to determine resolution. Numerical simulations suggest that the SRT algorithm produces fast and accurate reconstructions at realistic noise levels. The contrast is over 95% in all phantoms examined and is independent of noise level.


Introduction
In classical tomography, one needs to reconstruct a two-dimensional (2D) function f (x 1 , x 2 ) from the following projections: where ρ = x 2 cos θ − x 1 sin θ and τ = x 2 sin θ + x 1 cos θ. The above line integral is known as the 2D Radon transform. Its inversion was first constructed by Radon [1], and it is equivalent to the filtered backprojection (FBP) reconstruction formula [2,3,4]. Filtered backprojection (FBP) is the predominant analytic reconstruction method today. The main advantages of FBP are speed and simplicity. FBP assumes a simple Radon model where the data consist of line integrals along the radioactivity distribution, ignoring the randomness of the gamma-ray counting process. Therefore, in FBP it is difficult to incorporate complex physical phenomena such as attenuation and scatter. Noise issues are treated by selecting appropriate filtering parameters, such as the roll-off and cutoff frequencies of the reconstruction filter (usually at the expense of spatial resolution). Another disadvantage of FBP is the streak artifacts that are particular prominent near hot regions of the object.
In this study, we present an improved numerical implementation of an analytic, twodimensional, reconstruction method called SRT (Spline Reconstruction Technique), which was presented earlier in the literature [5,6]. This technique involves the numerical evaluation of the Hilbert transform of the sinogram via an approximation in terms of cubic splines. In this new version we have corrected singularity issues and optimized reconstruction time by employing mathematical symmetries and imposing sinogram thresholding to restrict reconstruction within object pixels. We have implemented SRT in the Open Source software library called STIR (Software for Tomographic Image Reconstruction), and have evaluated this method using simulated data from a clinical PET system.
and denotes the principal value integral.
For the numerical calculation of the Hilbert transform of f (ρ, θ) we assume that f (ρ, θ) has support in the interval −1 ≤ ρ ≤ 1 with f (−1, θ) = f (1, θ) = 0, and that f (ρ, θ) is given for every θ at the n points {ρ i } n 1 . In the interval ρ i ≤ ρ ≤ ρ i+1 , we approximate f (ρ, θ) by cubic splines: where . Note that f i denotes the value of f at ρ i , and f ′′ i denotes the value of the second derivative of f (ρ, θ) with respect to ρ evaluated at ρ i . By substituting Eq. (4) in Eq. (3) we can derive the following expression for the partial derivative of F (ρ, θ) with respect to ρ: where −1 ≤ ρ ≤ 1 and 0 ≤ θ < 2π. C(θ) and {D i (ρ, θ)} n 1 are defined as follows: The quantities can be computed in terms of f i (θ) by solving the following n linear equations: We note that the points {ρ i+1 } n−1 i=1 are removable logarithmic singularities. This is a direct consequence of Eqs. (7a).
Therefore, the Inverse Radon transform of a function f (ρ, θ), using our splines approach, can be written as: where C(θ) and D i (ρ, θ) are given by Eqs. (6a) and (6b), respectively. We note that in the construction of the so-called 'natural' splines, one requires continuity of the first derivatives (the set of equations in Eq. (7a)), as well as the conditions f ′′ 1 = f ′′ n = 0. The former requirement implies that there cannot be logarithmic singularities at the interior points ρ = ρ i , i = 2, . . . , n − 1. In order to eliminate the logarithmic singularities at the end points ρ 1 = −1 and ρ n = 1, we impose the set of equations in (7b) (instead of f ′′ 1 = f ′′ n = 0). In this way we construct a set of splines 'custom made' for the evaluation of the Hilbert transform.
In order to numerically evaluate Eq. (8) for f (x 1 , x 2 ) from the given sinogram f (ρ, θ) and the known detector locations ρ and projection angles θ, we first evaluate the second derivatives of the sinogram f ′′ (ρ i , θ) by solving the system of linear equations in Eqs. (7a) and (7b). Then, we calculate the term C(θ) for all θ's and ρ i 's. Finally, for any x 1 , x 2 , we compute ρ for every θ, and then f (x 1 , x 2 ) using Eq. (8).
The reconstruction time of this algorithm can be reduced by employing object specific information that is 'hidden' in the sinogram. In this respect we consider the important case that the boundary of the object is convex. In this case, a pixel which is outside the boundary spanned by an object and hence has zero value, can be singled out from the sinogram by first identifying the detector locations for all angles θ that receive contribution from this pixel; then, for every (x 1 , x 2 ), if there is even one θ such that f (ρ, θ) = 0, it follows that f (x 1 , x 2 ) must be zero.
Using the above condition we can restrict the reconstruction process only to pixels within the object boundary and exclude all zero pixels outside the object. In this way, in addition to improving considerably the reconstruction time, we can also obtain a 'clean' reconstruction without any streak artifacts outside the object. Note that for real data, the condition f (ρ, θ) = 0 must be replaced by f (ρ, θ) ≤ threshold, since in the presence of system noise, pixels outside the object's boundary in the sinogram can have values greater than zero. The threshold can be determined manually by examining the sinogram values outside the object boundary, or automatically using thresholding selecting techniques.
We have implemented the above algorithm in STIR [7], which is an object-oriented library using C++. For this purpose, we have utilized STIR's built-in classes and also have created a new class to accommodate our algorithm. The speed of the algorithm was optimized by using certain mathematical symmetries which exist in the standard case of constant detector spacing, which is indeed the case for the Discovery ST PET scanner employed in this study.

Simulation studies
We have used STIR to simulate the GE Discovery ST PET scanner [8]. In order to evaluate the performance of the SRT algorithm we have employed an image quality (IQ) phantom, a Jaszczak phantom and a slice of the digital 3D Hoffman phantom. The IQ phantom simulates the human torso and it consists of two circular cold regions (with diameter of 38 mm and 32 mm) and four circular hot regions (with diameters of 25 mm, 19 mm, 15 mm, and 12 mm) inside a larger warm region that simulates background. The radioactivity ratio between hot regions and the surrounding background is 4:1 for the three hot regions. The Jaszczak phantom is separated into six sections. Each section has cold lesions of different size uniformly arranged to form an equilateral triangle (with circle diameters of 27 mm, 18 mm, 16 mm, 13 mm, 11 mm, 9 mm) inside a hot region 30.67 cm in diameter. The center-to-center distance between lesions of equal diameter is twice their diameter. The Hoffman phantom simulates the radioactivity distribution of a cerebral PET study and it consists of three distinct radioactive regions: Gray Matter (GM), White Matter (WM), and Cerebrovascular Fluid (CSF).
After placing the three phantoms in the center of the scanner, 2D projection data were generated in STIR using a ray tracing technique with 10 rays per detector. Scatter and attenuation were not modeled. These sinograms provide the noiseless PET measurements. For each noiseless sinogram, 20 Poisson noise realizations have been generated at 3 different levels (NL1-NL3), where NL3 corresponds to the highest noise level. No filtering or smoothing has been applied to the SRT reconstructed images post reconstruction.
The contrast for the hot and cold regions were calculated according to [8]. For the IQ phantom, the contrast for the three smallest hot lesions was determined in each noise level by drawing an appropriate size circular ROI to determine the average counts measured in hot lesion and background, respectively. These values were averaged over all noise realizations. Similar procedure was followed for the Jaszczak phantom, where the contrast for the 18 mm, 16 mm and 13 mm cold lesions was determined by calculating the mean activity in all lesions of the section. Furthermore, a line profile was drawn through the smallest lesions of the Jaszczak phantom at NL3, in order to determine the capability of SRT to distinguish two cold lesions placed close to each other. The line profile has been normalized to its maximum value.

Results
The reconstruction time for SRT was 2.1 sec per sinogram, executed on a PC with Intel R ⃝ Core TM i7-920 Processor. We note that no parallel programming or other accelerating techniques have been employed.
Visual comparisons of the reconstructed images with no noise, as well as with noise are shown in Fig. 1. The noisy images presented are representative reconstructions of one realization at the specific noise level. Note that the SRT algorithm can generate negative values in the reconstructed images. In all images presented, the all-black color corresponds to zero values, whereas the white color represents the maximum value of the distribution. No streak artifact are present outside the object boundary.
The contrast for the three smallest hot spheres of the IQ phantom (19 mm, 15 mm, and 12 mm) as a function of noise level, is presented in Fig. 2A. Furthermore, the contrast for three cold sections for the Jaszczak phantom (with diameters of 16 mm, 13 mm, 11 mm) is illustrated in Fig. 2B. The SRT algorithm exhibits high contrast (over 95%) in all cold and hot lesions independently of noise level. Especially, SRT delivers almost 100% contrast for the cold regions.

Conclusion
In this work, we have presented the mathematical formulation and numerical implementation of a new, analytic, 2D, image reconstruction method for parallel beam geometry. Simulated images  have been reconstructed for various noise levels and have been evaluated in terms of cold and hot lesion contrast. Overall, the SRT algorithm provided high contrast images with no streak artifacts. Furthermore, unlike other analytic reconstruction algorithms, the reconstruction time of SRT is comparable with that of FBP.