Fast and Accurate Collimator-Detector Response Compensation in High-Energy SPECT Imaging with 1D Convolutions and Rotations
Lucas Polson Pedro Esquinas Sara Kurkowska Chenguang Li Peyman Sheikhzadeh Mehrshad Abbassi Saeed Farzanehfar Seyyede Mirabedian Carlos Uribe Arman Rahmim

 
Fast SPECT CDR Modeling]Fast and Accurate Collimator-Detector Response Compensation in High-Energy SPECT Imaging with 1D Convolutions and Rotations



^1 Department of Physics & Astronomy, University of British Columbia, Vancouver, Canada
^2 Department of Integrative Oncology, BC Cancer Research Institute, Vancouver Canada
^3 Molecular Imaging and Therapy Department, BC Cancer, Vancouver Canada
^4 Department of Nuclear Medicine, Pomeranian Medical University, Szczecin, Poland
^5 Nuclear Medicine Department, IKHC, Faculty of Medicine, Tehran University of Medical Science, Tehran, Iran
^6 Department of Radiology, University of British Columbia,Vancouver, Canada

lukepolson@outlook.com

§ ABSTRACT 
Objective: Modeling of the collimator-detector response (CDR) in SPECT reconstruction enables improved resolution and more accurate quantitation, especially for higher energy imaging (e.g. ^177Lu and ^225Ac). Such modeling, however, can pose a significant computational bottleneck when there are substantial components of septal penetration and scatter in the acquired data, since a direct convolution-based approach requires large 2D kernels. The present work presents an alternative method for fast and accurate CDR compensation using a linear operator built from 1D convolutions and rotations (1D-R). To enable open-source development and use of these models in image reconstruction, we release a SPECTPSFToolbox repository for the PyTomography project on GitHub. Approach: A 1D-R CDR model was formulated, and subsequently fit to Monte Carlo 440 keV point source data representative of emissions in ^225Ac imaging. Computation times of (i) the proposed 1D-R model and (ii) a traditional model that uses 2 dimensional convolutions (2D) were compared for typical SPECT matrix sizes. Both CDR modeling techniques were then used to reconstruct ^225Ac phantom and patient data, and were compared by quantifying total counts in hot regions of interest (ROIs) and activity contrast between hot ROIs and background regions. Results: The 1D-R and 2D CDR models were created using the SPECTPSFToolbox. For typical matrix sizes in SPECT reconstruction, application of the 1D-R model provides a two-fold computational speed-up over the 2D model. Only small differences between the 1D-R and 2D models (order of 1% ) were obtained for count and contrast quantification in select ROIs. Significance: A technique for CDR modeling in SPECT was proposed that (i) significantly speeds up reconstruction times, and (ii) yields nearly identical reconstructions to traditional 2D convolution based CDR techniques. The released toolbox  will permit open-source development of similar models for different isotopes and collimators.



















§ INTRODUCTION


Single photon emission computed tomography (SPECT) imaging is an in vivo modality of significant value in different clinical applications <cit.>. In particular, quantitative SPECT imaging holds great potential and value in the field of theranostics, wherein image quantification and dosimetry can be performed to improve radiopharmaceutical therapies (RPTs). Routine SPECT-based dosimetry could unveil the relationships between tumor-response and healthy organ complications with absorbed dose, therefore enabling personalized treatments that maximize tumor dose while minimizing toxicity to organs at risk <cit.>. As a recent and well-known example, the success of ^177Lu-PSMA-617 based RPTs in the VISION <cit.> and TheraP <cit.> randomized control trials has lead to expanded research and clinical translation efforts for wide-scale deployment, as well as pursuit of absorbed dose estimation as a tool towards personalized treatments <cit.>. 

SPECT-based dosimetry relies on accurate quantification of activity distribution from SPECT images, which is strongly dependent on the system model used in image reconstruction. The system model is a linear operator that predicts the expectation of the acquired data given a 3D isotope distribution, accounting for phenomena such as photon attenuation in the patient. An important aspect of the model is collimator detector response (CDR) modeling, which estimates image blurring caused by the collimator and the scintillation crystals <cit.>. Computation of the CDR is typically the most computationally expensive operation in image reconstruction.

The CDR can be characterized by the detector response from a point source of activity; this is denoted the point spread function (PSF). The PSF can be decomposed into multiple components. The intrinsic response function (IRF) characterizes the inability of the scintillation crystals to precisely localize the point of interaction, and is sufficiently modeled using a Gaussian function. The collimator response for parallel hole collimators results from the inability of the collimator to accept only photons travelling perpendicular to the detector; it consists of three components: (i) the geometric response function (GRF) <cit.>, which describes photons that travel through the collimator holes without penetrating or interacting with the septa, (ii) the septal penetration response function (SPRF), which describes the contribution from photons that travel through the collimator without being attenuated, and (iii) the septal scatter response function (SSRF), which consists of photons that interacted and scattered within the collimator and were subsequently detected in the scintillator. 

Selection of collimator parameters is an important aspect of SPECT imaging <cit.>. As the collimator becomes thicker and the diameter of the collimator bores becomes narrower, the relative contribution from the SPRF and SSRF decreases due to the increased attenuation probability for photons not travelling perpendicular to the detector. In this case, (i) the detector resolution improves and (ii) the net point spread function (PSF) is dominated by the IRF + GRF and can be reasonably approximated using a 2D Gaussian function.  A trade-off of having a collimator with high septal thickness and small hole diameters, however, is a decrease in detector sensitivity; the corresponding implications for quantitative imaging is a decrease in precision or longer patient scan times. This trade-off must be independently considered for each isotope. For ^177Lu labeled radiopharmaceuticals, where the imaged photons are 208keV, a commercially labeled “medium energy” collimator configuration (i) yields a reasonable count rate and (ii) adequately minimizes the SPRF and SSRF components. In this situation, a 2D Gaussian PSF model is sufficient for image reconstruction <cit.>.

Many commercially available reconstruction software only offer Gaussian PSF modeling, and thus implicitly assume the SPRF and SSRF are negligible. There are multiple advantages to this. Firstly, the computational advantage of 2D Gaussian PSF modeling is that 2D Gaussian convolution is separable into use of two perpendicular 1D Gaussian convolutions, which are less computationally expensive to implement. Secondly, the Gaussian PSF model used to model the GRF+IRF can be obtained analytically for any photon emission energy and standard collimator shapes, and thus does not require lookup tables. However, when the PSF has significant contributions from the SPRF and SSRF, 2D Gaussian PSF modeling fails to capture all the features of the PSF. This is typically an issue with radioisotopes that emit high photon energies, such as α-emitters like ^225Ac.

Use of α-emitters in radiopharmaceutical therapies presents a major and exciting frontier, due to the high linear energy transfer (LET) associated with α emissions <cit.>. In preclinical studies, simultaneous treatment with ^177Lu-PSMA-617 and ^225Ac-PSMA-617 compared to ^177Lu-PSMA-617 alone resulted in significantly reduced tumor growth <cit.>. Kratochwil et al. <cit.> applied ^225Ac-PSMA-617 to patients with metastatic castration resistant prostate cancer (mCRPC) who previously exhausted ^177Lu-PSMA-617 treatment. Certain patients achieved a full response, with prostate specific antigen (PSA) levels in one patient decreasing from 419 ng/mL to below 0.1 ng/mL. At the time of writing, there are ongoing studies of ^225Ac based radiopharmaceuticals looking at dose escalation <cit.>, fractionation <cit.>, and safety and efficacy <cit.>. A recent meta analysis <cit.> elaborates further on these clinical trials, discussing toxicities and other challenges with these treatments. Quite recently, targeted α therapy with ^213Bi has been shown to reduce amyloid plaque concentration in male mice <cit.>, eluding to a potential treatment option for Alzheimer disease.

Throughout the decay chain of ^225Ac, the daughters ^213Bi and ^221Fr emit photons detectable within a SPECT system of 440keV and 218keV respectively. Unfortunately, even with the  commercially available collimators designed for high-energy photons, there are still significant SPRF and SSRF components present in the ^213Bi 440keV peak; a sample PSF is shown in Figure <ref>. This phenomena similarly occurs in imaging of photons with energies of 511 keV (positron emitter) <cit.> and 364 keV (^131I) <cit.>. As a consequence, the PSF can no longer be modeled using a 2D Gaussian and is thus no longer separable into two perpendicular 1D components. Since this is the bottleneck of system modeling, image reconstruction takes significantly longer. 

The reduction of reconstruction times in medical imaging remains an important research topic. Tsai et al. <cit.> recently showed that the limited-memory Broyden-Feltcher-Goldfarb-Shannon algorithm with box constraints and a diagonal preconditioner (L-BFGS-B-PC) was able to converge several times faster than the one step late expectation maximum (OSL-EM) algorithm in Positron Emission Tomography (PET) reconstruction because less projection operations were required. As another example, in Chun et al. <cit.>, long reconstruction times for PSF modeling of ^131I were partially remedied by using fast fourier transform (FFT) based 2D convolutions. By contrast, in the present work, we seek to reduce reconstruction times by using a PSF model that incorporates (i) 1D convolutions and (ii) rotations. It will be shown that this model significantly reduces computational time compared to standard and FFT based 2D convolutions, and consequently significantly reduces the time required for SPECT image reconstruction of ^225Ac. The method is evaluated on phantom data and patient data, and we show that the proposed model yields identical results to Monte Carlo based 2D convolution techniques. All models considered in this paper are implemented using the GPU-accelerated functionality of PyTorch.

Alongside this paper, we release the SPECTPSFToolbox: an open-source GitHub repository which forms a new component of the PyTomography <cit.> project. The toolbox contains functionality for developing and fitting arbitrary PSF models to arbitrary point-source data. The saved models can then be loaded in our in-house initiated and community developed library PyTomography for SPECT reconstruction. As such, the techniques used in this paper can also be applied to other isotope / collimator configurations that are of interest in the nuclear medicine community. To encourage community use, we have released nine tutorials demonstrating how to use the toolbox, and how to integrate the models in PyTomography. The link to the PyTomography project (which includes the SPECTPSFToolbox) is <https://github.com/PyTomography>



§ MATERIALS AND METHODS


In Section 2.1, the mathematical formalism of the PSF model is outlined, and the reconstruction protocol for all acquired data in subsequent sections is established. Section 2.2 describes how the model is fit to SIMIND <cit.> Monte Carlo ^225Ac point source data acquisitions simulated at various source-detector distances. The computational time of the model is benchmarked and compared to 2D PSF modeling and Gaussian PSF modeling. In Section 2.3, the model is then used for reconstruction of (i) an ^225Ac phantom consisting of spherical inserts and (ii) a patient receiving ^225Ac-PSMA-617 treatment for metastatic prostate cancer. All computation was performed using a Microsoft Azure virtual machine (Standard NC6s v3) with a 6 CPUs (Intel(R) Xeon(R) CPU E5-2690 v4 @ 2.60GHz), 112 GB of RAM, and a Tesla V100 GPU.



 §.§ Theory




  §.§.§ PSF Modeling


In this work, the notation f(x,y;d) is used for 3D objects: d denotes the distance between a plane parallel to the detector and the detector, and (x,y) denote the position on the plane. It is assumed that x, y, and d are discrete and thus f consists of voxels. The following notation is used to represent convolution operator K:


    K^(x;d;b)f ≡∑_x' f(x',y;d) k(x-x';d;b)
    
        K^(y;d;b)f  ≡∑_y' f(x,y';d) k(y-y';d)
    
        K^(x,y;d;b)f ≡∑_x',y' f(x',y';d;b) k(x-x',y-y';d;b)

where b are additional parameters the kernel k depends on, such as collimator septal thickness L_b, hole diameter w_b, and the linear attenuation coefficient of the collimator material μ_b.

Assuming a linear shift invariant (LSI) PSF, the SPECT system matrix estimates the projection g_ϕ at angle ϕ as


    g_ϕ(x,y) = ∑_d∑_x',y' k_PSF(x-x',y-y',d) p_att(x,y,x',y',d) f(x,y,d)

where (x',y') is the source position on a plane parallel to the detector in 3D space, psf(...,d) is a kernel that yields the point spread function (PSF) at a distance d from the detector, and p_att(x,y,x',y',d) is the probability that photons traveling from (x,y,d) to detector coordinate (x', y') are not attenuated. Under the assumption that the attenuation probabilities vary little across the PSF, it follows that p_att(x,y,x',y',d) ≈ p_att(x,y,x,y,d) and it can simply be expressed as p_att(x,y,d). Defining f'≡ p_attf as the “attenuation-adjusted” image, Equation <ref> can then be rewritten in operator form (Equation <ref>) as


    g_ϕ(x,y) = ∑_d K_PSF^(x,y;d;b) f'

Since this convolution operation is typically the bottleneck of SPECT system matrix modeling and image reconstruction, it is of interest to look for techniques to reduce the computation time. Under conditions of no septal penetration and scatter, the CDR is dominated by the GRF: K_PSF^(x,y;d;b) can then be sufficiently approximated using a 2D Gaussian convolution K_G^(x,y;d;b) where the kernel is given by


    k_G(x,y;d;L_b,w_b,μ_b) = 1/2 πσ(d,L_b,w_b,μ_b)^2exp(-x^2+y^2/2σ(d,L_b,w_b,μ_b)^2)


    σ(d,L_b,w_b,μ_b) =  1/2√(2log(2))(w_b/L_b-2/μ_b· d + w_b )


Convolution with a Gaussian function has computational advantages since it can be decomposed into successive application of two perpendicular 1D kernels via K_G^(x,y;d;b)=K_G^(x;d;b)K_G^(y;d;b); 1D convolution is significantly more computationally efficient than 2D. Unfortunately, the decomposition of a 2D kernel into two 1D kernels is not mathematically possible when the CDR contains significant SPRF and SSRF components. 2D convolution, however, is not the only way to implement Equation <ref>. Owing to the discrete rotational symmetries and features of the anisotropic PSFs obtained with typical SPECT collimators, we propose the following “1D rotation” model (abbreviated as “1D-R”)  for the PSF operator:


    M_1D-R^(x,y;d;b)≡(∑_θ∈Θ_tℛ^-1_θ K^(x;d;b)_t ℛ_θ + K_B^(x;d;b) K_B^(y;d;b) + 1) K^(x;d;b)_G K^(y;d;b)_G

where ℛ_θ is a rotation operator that implements rotation about an axis in the d direction by angle θ, and Θ_t correspond to the angle of the septal penetration tails (equal to {0,π/3,2π/3} for ^225Ac). The derived linear operator only makes use of 1D convolutions and rotations, since they often require less computational time on GPU compared to 2D convolutions. The kernels, modeling the components of the CDR (as shown in Figure <ref>), are selected as follows:


    
  * Gaussian kernel K^(x;d;b)_G:
        
    k_G(x;d;b) = A_G(d,b) ·exp(-x^2 / 2σ_G(d,b)^2 )

        
    A_G(d,b) = b_0e^-b_1d + b_2e^-b_3 d

        
    σ_G(d,b) = b_4 + b_5(√(d^2+b_6^2) - |b_6|)

    This kernel is used to build geometric component of the PSF.
        
    
  * Tail kernel K^(x;d;b)_t:
        
    k_t(x;d;b) = A_t(d,b) f_t(x/σ_t(d,b))

        
    A_t(d,b) = b_0e^-b_1d + b_2e^-b_3 d

        
    σ_t(d,b) =1+b_4(√((d-d_min)^2+b_5^2) - |b_5| )

    where d_min is the source-detector distance used in the PSF fit. f_t(x) is a discrete array of numbers that is linearly interpolated between its fixed points; these points are also considered to be part of the hyperparameters b. This kernel is used to build the septal penetration component of the PSF.
    
    
  * Isotropic background kernel K^(x;d;b)_B:
        
    k_B(x;d;b) = A_B(d,b) f_B(x/σ_B(d,b))

        
    A_B(d,b) = b_0e^-b_1d + b_2e^-b_3 d

        
    σ_B(d,b) =1+b_4(√((d-d_min)^2+b_5^2) - |b_5| )

        
    f_B(x) = e^-|x|

        This kernel is used to build the septal scatter component of the PSF.




All three model components also encapsulate the small blurring contribution from the IRF. The hyperparameters b are separate for each of the three different components.

In subsequent sections, the proposed 1D-R model is evaluated against a Monte Carlo 2D model where kernel data is obtained using SIMIND <cit.>; point source projection data is normalized at a large number of source-detector distances. The corresponding PSF operator, abbreviated as “2D”, is defined as 


    M_2D^(x,y;d;b)≡ K_true^x,y;d

where the corresponding kernel k_true^x,y;d corresponds to normalized projection data at each source-detector distance d.

The publicly SPECTPSFToolbox repository is structured to facilitate the customizability of these models. For the 1D-R model, components are implemented via separate class instances that are subsequently added and multiplied together to obtain the final form of Equation <ref>. The functional form and hyperparameters of the amplitude/scaling, such as that of Equation <ref>, are left arbitrary for the user. The 2D model (Equation <ref>) is implemented via a class instance that receives a stack of PSF kernels with associated source-detector distances: when used in image reconstruction, the model automatically chooses the PSF kernel closest to each source-detector distance in the reconstruction problem.



  §.§.§ Image Reconstruction


All acquired image data in subsequent sections are reconstructed using PyTomography with the maximum likelihood expectation maximization (MLEM) algorithm; this was chosen over the ordered subset expectation maximization algorithm (OSEM) because of the low count statistics for each projection. System matrices employed attenuation correction using attenuation maps derived from acquired computed tomography (CT) images. As described by Zekun et al. <cit.>, SPECT with α-RPTs often have significant stray radiation-related noise component due to the proportionally low count rate. The imaging system equation is thus given by


    y̅ = Hx + Ψ̅ + ŝ


where y̅ is the expectation of the acquired count data y, H is the system matrix, x is the estimated 3D count distribution, Ψ̅ is the expectation of the stray radiation-related noise, and ŝ is the scatter estimate. Ψ̅ is assumed to be the same for all bins in a given energy window and proportional to the acquisition time, and can be measured by taking a blank scan using the same energy windows as the patient acquisition, averaging the number of counts in each valid bin, and scaling by the acquisition time used in the clinical protocol. Since the dual energy window (DEW) technique is used for scatter estimation, the stray radiation in that window also must be accounted for: the scatter estimation is thus given by


    ŝ = ( w_p/w_l) (S_G[g_l] - Ψ̅_̅l̅)


where w_p is the width of the primary window, w_l is the width of the lower window, S_G is a Gaussian smoothing kernel, g_l is the acquired count data in the lower energy window, and Ψ_l is the expectation of the mean stray-radiation noise in the lower energy window. The smoothing kernel is used here due to the high noise present. One of the dangers of Equation <ref> is that it is capable of yielding negative values in Equation <ref>; in practice this rarely occurred, and all negative values were set to zero.

When used in image reconstruction, the size of the 2D PSF kernels was always set to N-1 × N-1 where N is the number of voxels along the largest direction in the reconstructed image. For the 1D-R model, the kernels corresponding to K_G and K_B were of size N-1, while the kernel corresponding to K_t was of size ⌈√(2) N ⌉ + n to account for the diagonal, where n is of integer value 0 or 1 to make the kernel size odd. It should be noted that for PSFs much larger than the dimensions of the image, the size of the kernels would need to cover a 2N-1 × 2N-1 area to account for contributions from voxels on one edge of the image to detector elements on the opposite edge; in practice, the PSFs were not that large and the dimensions chosen were sufficient.



 §.§ Validation and Timing


A 440 keV point source was simulated using the SIMIND Monte Carlo program <cit.> at 1100 positions that linearly varied from 0 cm to 58.44 cm. The detector pixel size was 0.24 cm× 0.24 cm with 255×255 pixels. The simulated collimator corresponded to a Siemens high energy configuration with a hexagonal shape, a hole width of 0.4cm and a hole length of 5.97cm. The intrinsic resolution of the detector was also included in the simulation and was assumed to be 0.38cm at 140 keV, representative of a Symbia system with 3/8” crystal length.  The PSFs at all 1100 distances were used to generate K_true^x,y;d in Equation <ref> during reconstruction by selecting the PSF closest to the actual source-detector distance.

Twelve of these PSFs at source detector distances of 1cm and every 5cm from 5cm to 55cm were used as data to fit the 1D-R model. Fitting was performed in three steps: (i) the Gaussian parameters were optimized using the Adaptive Moment Estimation (ADAM) algorithm for 1 · 10^4 iterations and a learning rate of 10^-2, (ii) all the parameters were simultaneously optimized using ADAM for 1.5 · 10^4 iterations and a learning rate of 10^-3.

The computational time for the 1D-R model (Equation <ref>) and 2D model (Equation <ref>) were bench-marked on CPU and GPU for by applying the operators to matrices of four different sizes: 64^3, 128^3, 196^3, and 256^3. Each matrix of size N^3 was filled with random uniform numbers between 0 and 1. Each experiment used conditions similar to SPECT imaging, where the d axis varied from 0cm to 50cm, and the x and y axes varied from -30.74cm to 30.74cm. Use of standard convolution and fast fourier transform (FFT) based convolution were compared for each method. 



 §.§ SPECT Studies




  §.§.§ Phantom Study


Application of (i) standard Gaussian (ii) 1D-R (iii) 2D PSF modeling was evaluated for GPU reconstruction of acquired ^225Ac data. A cylindrical phantom with spheres of diameter 60mm, 28mm, and 22mm with an initial sphere activity concentration of 1.37 kBq/mL was filled at a 10:1 source to background ratio. 34 SPECT acquisitions of the phantom were taken in sequence on a Symbia T2 SPECT/CT system (Siemens Healthineers, USA) with the following settings: 128 × 128 pixels at 4.82 mm× 4.82 mm resolution, 96 projection angles, high energy collimators, and 60 s acquisition time per projection. For reconstruction of the 440keV emission from the ^213Bi daughter, two energy windows were configured for acquisition: (i) photopeak window centered at 440 keV with a width of 20% and (ii) a lower scatter window used for DEW scatter correction centered at 374keV with a width of 12.5%. A blank scan with identical parameters was acquired to obtain the mean-stray radiation noise in each energy window. Two noise levels of the data were considered:


    
  * One of the 34 scans was used for data reconstruction. The number of acquired counts in this scenario corresponds to an expected clinical ^225Ac scan of a patient undergoing ^225Ac-PSMA-617 therapy with injection of 8MBq <cit.> and is scanned for 2.5 min per projection at a time point between 0-72 hr post injection.
    
  * The counts from all 34 scans were summed together and used for image reconstruction. In this case, no scatter blurring was used in Equation <ref> since the noise level was low. Because the detectors had an approximately equal radial path around the phantom for each scan, this corresponds to a single scan with a high count rate.


Although 96 projection angles were acquired, only 32 were used since this was similar to the corresponding ^225Ac patient acquisition in the next section. All images were reconstructed with MLEM for up to 100 iterations.



  §.§.§ Patient Study


A patient receiving ^225Ac-PSMA-617 therapy with an injected activity of 8 MBq was imaged 20.5 hours post injection (first cycle). SPECT data were acquired on a Discovery 670 Pro SPECT/CT (GE Healthcare, USA) with a high energy general purpose collimator. 30 projections (15 per head) were acquired for 150 s using an identical photopeak window to the phantom acquisition, and 358.15 keV - 395.85 keV lower window. Since this GE scanner has a different CDR than the Siemens scanner used in prior examples, PSF data were regenerated in SIMIND using a high energy general purpose GE collimator and another 1D-R model was fit as before. In GE scanners, the hexagonal bores in the collimator are rotated by 90^∘ compared to the Siemens collimator, so the angles in Θ_t were adjusted to compensate for this. The acquired data were reconstructed using (i) standard Gaussian (ii) 2D and (iii) 1D-R PSF modeling using MLEM for up to 100 iterations. Reconstructed images were filtered using a Gaussian filter with a full width half max (FWHM) of 3cm. Two high uptake bone lesions were segmented on the 100th iteration reconstruction using 3D Slicer <cit.>. Three lesion ROIs were segmented by a physician on a pre-therapeutical PET image using the PET Edge+ tool of MIM v7.2.1 (MIM Software Inc., USA). A spherical ROI of diameter 10cm was placed in a low uptake region in the center of the patient so that the contrast, defined as the mean uptake ratio between the lesion ROIs and background ROI, could be obtained. The mean number of counts and contrast for each lesion ROI were then evaluated for each iteration of MLEM; before the statistics were computed, the 3cm FWHM Gaussian filter was applied to the image.



§ RESULTS


Figure <ref> shows the Monte Carlo simulated PSF data compared to the fit obtained via Equation <ref>. The curve fit is a reasonable approximation to the PSF at most distances, but is unable to capture features at larger distances, such as the intensity of the tails at 50 cm. The model also consistently under predicts the intensity at the center of the PSF, and this effect is more prevalent at large distances.   


Timing benchmarks for PSF modeling using the 1D-R, 1D-R (FFT), 2D, and 2D (FFT) are shown in Figure <ref>. CPU implementation yields no benefits with the proposed 1D-R model, and performs fastest with FFT-based 2D convolutions and slowest using regular 2D convolutions. The GPU implementation is faster than the CPU implementation for all methods, and yields computational benefits when the proposed 1D-R model is used. With a matrix size of 128^3, use of the 1D-R method is over three times faster than 2D (FFT) PSF modeling. For small matrix sizes (64^3) there is no computational speed-up using the 1D-R method, and for large matrix sizes (256^3) the relative time difference between the 1D-R and 2D (FFT) begins to decrease, approaching only a two times speed advantage.



PyTomography reconstructed ^225Ac phantom images are shown in Figure <ref>. The time required for 100 iterations of MLEM for the low count data was 110s (Gaussian PSF), 146s (1D-R PSF) and 377s (2D PSF). Reconstructions using the 1D-R and 2D methods are almost indistinguishable qualitatively. All 34 low count acquisitions were then reconstructed and the 1D-R and 2D methods were quantitatively compared; the differences between the 2D and 1D-R were (1.25 ± 1.25)%, (0.60 ± 0.69)%, and (1.90 ± 0.18)% from smallest to largest sphere respectively. For comparison, the variability in counts between separate noise realizations for the 2D model were 54.4%, 33.0%, and 9.8% from smallest to largest sphere, respectively.



PyTomography reconstructed ^225Ac-PSMA-617 patient images are shown in Figure <ref>. The time required for 100 iterations of MLEM was 38.5 s (Gaussian PSF), 165.6s (1D-R PSF) and 351.5s (2D PSF). Qualitatively, the Gaussian PSF yields less counts in the uptake ROIs, and also greater background counts, for example, in the vertebrae. The count differences between the 1D-R and 2D methods after 100 iterations were 1.27% (lesion 1), 0.48% (lesion 2) and 2.03 % (lesion 3). These differences are comparable to those observed in the phantom study, and were significantly less than the differences between the Gaussian and 2D methods of -42.3 % (lesion 1), -29.1 % (lesion 2) and -44.3 % (lesion 3). The differences in lesion-background contrast between the 1D-R and 2D methods after 100 iterations were -0.18 % (lesion 1), -0.96 % (lesion 2) and 0.56 % (lesion 3). These were smaller in magnitude than the contrast differences between the Gaussian and 2D methods of -51.3% (lesion 1) and -23.3 % (lesion 2) and -56.7 % (lesion 3). Based on the line plots in Figure <ref>, the Gaussian PSF yields convergence with less iterations; this observation is similar to how an absence CDR modeling leads to faster convergence in traditional low and medium energy SPECT.  In general, the smaller the spatial extent of the PSF, the faster convergence is attained in MLEM.





§ DISCUSSION


In this work, a 1D-R PSF model for high energy SPECT reconstruction was proposed, and demonstrated to be both efficacious and fast for use in high-energy reconstruction of ^225Ac data. The model was created using publicly shared SPECTPSFToolbox component of the PyTomography project; t is hoped that the shared toolbox and corresponding tutorials will permit the nuclear medicine community to develop custom PSF models for new isotope/collimator configurations as research in novel RPTs and imaging techniques continues to grow.

As shown in Figure <ref>, the proposed 1D-R only yields computational benefits when image reconstruction is performed on GPU. If reconstruction is performed on CPU, then PSF modeling should be implemented with FFT-based 2D convolutions. FFT-based 2D convolutions are still, however, three times faster on GPU than on CPU. Furthermore, if a GPU is used then the 1D-R model can be used for further computational advantages. The fastest implementation of PSF modeling on GPU (proposed 1D-R) is approximately one order of magnitude faster than the fastest implementation of PSF modeling on CPU (FFT-based 2D convolutions).

While the proposed 1D-R model is not able to perfectly capture the features of the PSF at all radial distances in Figure <ref>, it yields near identical reconstructions in Figures <ref> and <ref>, especially when compared to reconstruction using the Gaussian kernel. In the low count scenario of the phantom study, which approximately represents a clinical count rate, the percent difference in the three spheres between the 1D-R and 2D model were small compared to the variability between separate acquisitions. In the patient example, the variability of counts between the 2D and 1D-R models in the bone lesions was similar to the spheres in the phantom: this suggests that any differences in predicted counts between the 2D and 1D-R models in the patient example is also negligible compared to differences that would be observed between separate acquisitions. While the 1D-R model and 2D model yielded similar contrast, the Gaussian model had significantly reduced contrast with differences of -51.3 %, -23.3 %, and -56.7 % compared to the 2D model in the bone lesions. While scaling the magnitude of the image reconstructed using the Gaussian PSF could artificially increase the number of counts in each bone ROI, it would have no impact on these substantial differences in contrast.

In the phantom example, the radial path of the detector approached 7.2cm at its minimum and 27.2cm at its maximum. Since the radius of the cylinder was 10cm, this suggests that the PSF model was able to reasonably function between 0-37cm. Future work may seek to test the efficacy of the proposed 1D-R model for larger radial distances, as this may be required for larger patients.

While the functional form of Equations <ref>, <ref>, <ref> were selected because they produced a reasonable approximation of the PSF, it may be the case that they can be improved upon. Use of rotation operations with the component K_B^(x;d;b)K_B^(y;d;b) may permit more radial symmetry in the SSRF. Substitution of K_G^(x;d;b)K_G^(y;d;b) with a small 2D kernel that represents the true aperture function of the GRF may remedy the under prediction in the center of the PSFs at large distances. Close inspection of the SIMIND data at large distances in Figure <ref> reveals additional dim tails; these could be modeled by including another independent tail component in the model. A trade-off for adding additional features to the model, however, is an increase in computation time. The SPECTPSFToolbox Python library released along with this paper contains tutorials demonstrating how to create custom and fit parameterized PSF models. Users can obtain point source data using Monte Carlo programs (such as GATE <cit.> and SIMIND) or by using real scanners, and fit corresponding models to the acquired data. The models can then be imported to PyTomography for customized and scanner specific SPECT reconstruction. Users can independently evaluate the trade-off between adding model features, impact on reconstructed images, and increase in computation time.

As discussed in Section 2.1.2, the kernel sizes were fixed to N-1 × N-1 for all d. For small d, where the PSF is also smaller, a kernel of this size might be excessive and could unnecessarily increase the computational time required. Future research exploring the reduction of computational time might experiment with using a source-detector distance dependent kernel size N(d) for that matches the size of the PSF; this would further reduce the computational time required. 

Since reconstruction using the proposed 1D-R model more than halves reconstruction times compared to the 2D model, it may be preferable for reconstruction of ^225Ac data. Both the 1D-R and 2D models, however, are only applicable when Equation <ref> is used for SPECT system matrix modeling; this equation relies on the assumption that the attenuation probabilities vary little across the PSF. Due to the large spatial extent of the PSFs shown here, this assumption may be invalid, and Equation <ref> may instead be required for accurate reconstruction of ^225Ac data. This assumption largely depends on the density of the object being scanned. If the object had a region of high density which only rays on the outer edges of the PSF passed through, then the equal attenuation path assumption would be invalid. For patients and phantoms, where the density is roughly constant throughout the field of view, application of Equation <ref> may be permitted. Meanwhile, the study of different density configurations and the applicability of Equation <ref> remains to be investigated in future work.

In conclusion, a fast and accurate implementation of high energy CDR modeling in SPECT imaging that uses 1D convolutions and rotations was developed, and tested on ^225Ac reconstructions. The technique was implemented using the open source reconstruction library PyTomography, and was shown to speed up reconstruction times by more than a factor of two compared to conventional, 2D convolution based methods. Furthermore, the 1D-R method yielded near identical results to the 2D method, and the small differences were insignificant compared to the differences between reconstructions of separate noise realizations. The SPECTPSFToolbox python library was developed and publicly shared to permit others in the nuclear medicine community to develop custom isotope/collimator PSF models for use in the open-source reconstruction library PyTomography.



§ ACKNOWLEDGEMENTS


This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) CGS D Award 569711 and Discovery Grant RGPIN-2019-06467, as well as computational resources and services provided by Microsoft AI for Health.



§ REFERENCES


unsrt 






