|Home | About | Journals | Submit | Contact Us | Français|
Monte Carlo (MC) is a well-utilized tool for simulating photon transport in single photon emission computed tomography (SPECT) due to its ability to accurately model physical processes of photon transport. As a consequence of this accuracy, it suffers from a relatively low detection efficiency and long computation time. One technique used to improve the speed of MC modeling is the effective and well-established variance reduction technique (VRT) known as forced detection (FD). With this method, photons are followed as they traverse the object under study but are then forced to travel in the direction of the detector surface, whereby they are detected at a single detector location. Another method, called convolution-based forced detection (CFD), is based on the fundamental idea of FD with the exception that detected photons are detected at multiple detector locations and determined with a distance-dependent blurring kernel. In order to further increase the speed of MC, a method named multiple projection convolution-based forced detection (MP-CFD) is presented. Rather than forcing photons to hit a single detector, the MP-CFD method follows the photon transport through the object but then, at each scatter site, forces the photon to interact with a number of detectors at a variety of angles surrounding the object. This way, it is possible to simulate all the projection images of a SPECT simulation in parallel, rather than as independent projections. The result of this is vastly improved simulation time as much of the computation load of simulating photon transport through the object is done only once for all projection angles.
The results of the proposed MP-CFD method agrees well with the experimental data in measurements of point spread function (PSF), producing a correlation coefficient (r2) of 0.99 compared to experimental data. The speed of MP-CFD is shown to be about 60 times faster than a regular forced detection MC program with similar results.
Image degradation effects, such as photon scatter, attenuation, and septal penetration, can greatly affect the quality of single photon emission computed tomography (SPECT) images , . Because of the complexity of these effects, corrections for them are typically inexact and may, in fact, cause artifacts in the reconstructed images –. The use of more accurate modeling of the acquisition process in SPECT can result in reduced artifacts, improved detection ability, and more accurate quantitation of the reconstructed image . Monte Carlo (MC) simulation of photon transport is one of the well-established tools that has been used in SPECT image reconstruction due to its ability to accurately model photon transport , . However, the high computational demand of MC limits its use primarily to research. In standard MC simulation of SPECT imaging, it is common that only one out of 105 emitted photons is actually detected by the gamma camera . Therefore, in order for MC modeling to be routinely used for clinical practice, it is necessary to improve the simulation speed , . Several variance reduction techniques (VRT) have been applied in order to decrease the computation time as well as increasing the photon detection efficiency. One such technique is called forced detection (FD) , . In FD, after emission of a photon from an initial location, the photon may scatter at various sites in the object prior to escaping. At each scatter site, the photons can be forced in the direction of gamma camera and incident on the detector, as shown in Fig. 1(a). The probability of detecting those photons are modified by multiplying the initial photon weight at the scatter sites with the probability that photon travels through the forced direction and is detected. The forced path direction is sampled from a probability density function (PDF) obtained from the model of distance-dependent collimator response. This technique dramatically increases the simulation speed compared to standard MC, but it is common for FD simulations to still take several hours or longer to achieve a noise-free projection image.
As an improvement over FD, de Jong et al. have developed a flexible technique called convolution-based forced detection to speed up FD . The fundamental idea of CFD is that the photon flux emitted from a point source and incident on the detector can be described as a function of the distance dependent response of the detector. That is, the flux is spread over the detector with approximately a Gaussian blurring kernel. This kernel is based upon the distance between the point source and detector, various collimator parameters, and the intrinsic resolution of the detector. Therefore, in the CFD process, unlike sampling the path direction from the PDF when forcing the detection, the photons are forced to travel along the path perpendicular to the detector while convolving with a Gaussian blurring function. In this method, the photon history weights are stored in a three-dimensional (3-D) sub-projection map before the value of each voxel is convolved with the Gaussian kernel and multiplied by the attenuation on the forced photon path corresponding to each voxel. It has been seen that this code has increased the speed more than 80 times compared to standard FD MC. However, this method may store several photon weights in the one pixel although their corresponding photon energy may be different, which may result in roughly calculated attenuation by using the average linear attenuation coefficient (μ) for each pixel. In order to achieve more accurate CFD simulation, we recently investigated another CFD code , as shown in Fig. 1(b). In this code, we direct photons at each location to the detector before calculating the attenuation from every photon location to the detector and, furthermore, using the attenuation coefficients based upon the individual photon energy. The sub-projection map is no longer required in our code because the photons are followed at each location until they escape the object. This code is developed in SIMIND MC and can be easily applied into other MC codes. Due to the additional overhead of calculating the exact attenuation path for each photon location, previous tests have shown that our CFD is about 10 times faster than FD .
de Jong et al. further developed a rotation-based MC (RMC) method in order to achieve several projection images from a few base projections . For more flexible CFD application, in this paper, we introduce another SPECT-type (i.e., rotating detectors) acquisition by modeling multiple projection measurements simultaneously. This multiple projection (MP) sampling is combined with the our CFD method (MP-CFD) in SIMIND MC simulation code in order to further reduce simulation times.
The following section will first explain the MP-CFD method in detail and how it is incorporated into the SIMIND MC program. The accuracy of this method will then be evaluated by comparing the results of it with experimental data and comparing the speed of it with more conventional CFD methods in order to measure any performance improvements.
In this section, we introduce the model of distance-dependent collimator response and intrinsic resolution upon which the CFD method is based, followed by the discussion of photon transfer probability calculation and the model of MP-CFD.
In SPECT imaging, the overall detector resolution is a combination of the intrinsic camera resolution Ri and the specific collimator resolution Rc. The intrinsic resolution refers to how precisely the scintillator and associated localization electronics can position an event to a specific location. It is normally a property of the camera and cannot be changed. The collimator resolution, on the other hand, can be altered as it is dependent on the collimator length l, hole size d, source distance z over the collimator, and linear attenuation coefficient μ for the collimator material (e.g., lead) at a certain photon energy. The overall detector resolution can be described in terms of full width half maximum (FWHM) :
Furthermore, Rc can be written as
Normally, l − 2/μ is called the effective collimator hole length as it accounts for some septal penetration. It has been shown that the geometric response [point spread function (PSF)] of a collimator can be modeled as a Gaussian function
where (x0, y0) is the central location where photons are detected on the collimator face and where σx and σy are derived from the FWHM of the particular collimator.
The SIMIND MC program is a well-used program to simulate the photon transport for SPECT imaging . This program builds an accurate mathematical model of the SPECT detector system, incorporating different photon interaction types such as photoelectric effect, coherent scattering, Compton scattering, and pair production.
With SIMIND, photons are emitted according to the emission distribution map with an initial photon detection weight set to unity. The photon direction, path length and interaction type are sampled by a series of random numbers R. At every interaction site, the probability for the transportation of this photon is calculated repeatedly until it escapes the object, and the sampling of another photon emission starts.
Various techniques have been introduced in order to accelerate the speed of MC simulation. Collectively, these techniques are referred to as variance reduction techniques (VRT). One of these techniques is FD, in which the emitted photons are directed towards the detector, thus resulting in the detection of each photon. The exact path of the forced direction is sampled from a probability density function that is a function of the collimator and the photon-detector distance. The resultant photon event probability is then corrected by the probability that the photon travels directly to the collimator and is absorbed by the crystal and is convolved with a δ function at the camera location (x, y), as seen in Fig. 1(a).
The convolution-based forced detection (CFD) method is very similar to FD, with the exception that the photons are forced to travel in the path perpendicular to the gamma camera, with the detection probability being convolved with a Gaussian function determined by (3), as illustrated in Fig. 1(b). As discussed previously, the convolution work is performed on every photon transverse location, and the attenuation coefficient of each directed path is calculated using a specific μ value corresponding to the given photon energy E.
The calculation of photon probability P(α) of a photon reaching a certain detector when the detector is located at the angle of α is very important in CFD. The main photon interaction types in SPECT are Compton scatter and photoelectric interaction. During the transport of each photon, the probability that a photon travels to every interaction position is determined by the product of photon history weight (PHW), which starts at 1 at the initial emission location, the probability that the photon is not absorbed by the photoelectric interaction being (1 − Ppe (E)), and the probability of scattering due to the Compton effect based upon different scatter angle ϕ. The Klein–Nishina (KN) function, which is denoted as KN (ϕ), is introduced to calculate the probability of Compton effect. As a result, a new photon energy (E) and direction are given for the further photon transmission. At each scattered location, the attenuation of the CFD-forced photon travel path will also influence the photon history.
The probability of a photon being detected by a detector at angle can be given by
where i is the photon interaction number.
In both FD and CFD, photons are forced to travel towards the gamma camera to be detected. It is noted that the detection of these photons occurs at only a single projection angle for each simulated photon. When a SPECT acquisition is simulated, however, the projection data must be acquired from cameras at different projection angles. Thus, when using conventional FD and CFD, the simulation time increases by the number of projection angles simulated.
We now introduce the multiple projection (MP) simulation method, which will simulate projection images at different angles around the object simultaneously. An illustration of this method is shown in Fig. 2(a) and (b). Several detectors are modeled simultaneously with the number of detectors and their angular distributions being dependent upon the SPECT acquisition parameters. As the photon is emitted from the site of decay to a certain interaction point, it is directed to each of the detectors surrounding the object. The CFD method is implemented to every detector by forcing the photon to travel along the path directly to each of them, which is depicted in Fig. 2(b).
When the photons are directed to different cameras, the probabilities are different at each angle due to differing attenuation paths and the probability of Compton scattering. Thus, the overall detection probability of detecting a photon at detector angle θ can be written as
where P(α) denotes the photon event detection probability for detector at the angle of α, which acts as the base projection detector. P(θ) is the photon event detection probability for detector at angle of θ, and is the forced path (FP) length in the direction of θ angle detector. μ(E) is the linear attenuation coefficient corresponding to the photon energy E. The term gives the attenuation along the path perpendicular to the detectors at the angle of θ from the photon scatter locations to the corresponding detectors. In practice, we assume the probability due to KN function has been included into the attenuation terms. Therefore, (6) can be approximated as
The bias due to the scattered photons has been discussed in FD. The MP-CFD method did not produce extra bias when simulating multiple detectors around the object.
In this section, we present the validation of MP-CFD by comparing it to normal CFD results and experimental data. A GE Millennium VG dual head gamma camera is simulated using high energy general-purpose (HEGP) collimator and I-131 isotope. The FWHM is calculated for each point spread function (PSF). The correlation coefficient (r2) is introduced to measure the difference in the FWHMs for all the PSFs generated by MP-CFD and experimental measurements. The r2 can be calculated by
where xi is the experimental FWHM value, is the mean value of x, yi is MP-CFD FWHM value, and is the mean value of y.
In order to evaluate the accuracy of 2-D projection images obtained using MP-CFD, the linear and log center profiles of PSF produced by MP-CFD are compared with the experimental data. Furthermore, we have used the normalized mean square error (NMSE) to verify the accuracy of MP-CFD with independently obtained CFD projection images as the standard images. Thus
where Oi and Fi are the ith pixel values of the MP-CFD projection images and the CFD projection images, respectively.
The speed increase of MP-CFD has been verified by comparison with the standard CFD based upon the SPECT acquisition of the 64 × 64 × 64 NURBS-based Cardiac-Torso (NCAT) phantom. We have chosen to compare the MP-CFD method with CFD at this point as we have previously shown that the CFD method achieves very similar results to the standard FD but is approximately 10 times faster . Note that in order to compare the detection efficiency of MC with the experimental data, the unit of the projection images is set as counts per second per MBq (CPS/MBq). The number of photon counts detected per second would not be equal to the total radioactivities because of the photon transmission directions, attenuation, scatter, and other factors. The MC code will simulate the entire system, which is similar to the real detectors. The convolution kernel of CFD is related to the photon probability function, with the normalized sum. Therefore, CFD technique will not affect the detection efficiency.
To test the accuracy of the MP-CFD methods, the projection images of a distance dependent point source are obtained in order to compare with the experimental measurement. HEGP collimator has been used with 200 MBq I-131 point source located both in air and water. Fig. 3 shows the experimental setup with point sources located in the attenuation medium (non-radioactive water). The acquisition time for each frame is 10 min. The point source is moved along the centered axis of the phantom, varying the source detector distance. The normal CFD method is run twice for each source location in order to acquire simulated images for both detector heads, whereas MP-CFD code only needs to simulate once to get the two images of both heads. However, in order to validate the effect of base projection, MP-CFD is performed twice using head 1 and head 2 as base projections, respectively. Therefore, MP-CFD will generate four projection images for each point source.
The center profile through each projection image has been fit to a Gaussian function. The FWHM of the Gaussian function and the total counts of each projection image detected per second (CPS) given 1 MBq activity are recorded to compare with the experimental data.
Fig. 4 is the measurement of point sources in air away from head 1 at 2.54, 7.62, 17.78, 30.48, and 43.18 cm. The distance from them to head 2 are 53.34, 48.26, 38.1, 25.4, and 12.7 cm, respectively. The result of MP-CFD given in Fig. 4 is the measurement of both head 1 and head 2, using head 1 as the base projection.
Fig. 4(a) and (b) depict the very good agreement of normal CFD, MP-CFD generated images detected by both detector heads with the experimental data by the horizontal and vertical FWHM measurements. Fig. 4(c) and (d) shows that septal penetration (SP) will result in about 8 CPS more for 1 MBq activity compared with non SP simulation.
Fig. 5 shows a comparison of the horizontal and vertical direction FWHM measurements along with the relative sensitivity detected by these two heads when the point source is located in water. The point source distances from head 1 and head 2 are 15.24, 22.86, 27.84, and 33.02 cm and 40.64, 33.02, 27.94, and 22.86 cm, respectively. The subsequent depths in water corresponded to 5.30, 12.92, 18.00, and 23.08 cm and 18.20, 10.58, 5.50 cm, and 0.42 cm. Fig. 5(a) evaluates the accuracy of MP-CFD when using head 1 as the base detector, and Fig. 5(b) shows the evaluation of MP-CFD when using head 2 as the base detector. The first two columns of Fig. 5(a) and (b) give the evaluation of MP-CFD generated images on head 1, with the other two columns showing the evaluation of the head 2 simulation. The top rows of Fig. 5(a) and (b) compare the detector resolution in terms of FWHM. For easier analysis, the x-axis has been set as the distance between point source to the detectors. It is not surprising that MP-CFD matches CFD for base projection images, and they both agree with the experimental data. For the off-base projection images, the difference between MP-CFD and CFD simulations and the experimental data is slightly different, which may result from the approximated updating of (7). The bottom rows of Fig. 5(a) and (b) provide the comparison of the average photon counts detected by these two heads. The x-axis has been set as the depth of water from the point source to the detector in order to observe the effect of attenuation. It has been noted that our simulation results compare well to the detection sensitivity of head 2 but not head 1. Table I gives the r2 value of FWHM measurement of PSF’s in both air and water; the values of FWHM of PSFs in water is the minimum r2 value of the corresponding terms in Fig. 5. The result of MP-CFD agrees well with CFD and experimental data as r2 values for FWHM are all higher than 0.99.
In order to further verify the accuracy of the PSFs generated by MP-CFD, Fig. 6 presents the profiles of the two projection images detected by head 1 and head 2 using head 1 as the base projection, when the point source is placed 22.86 cm away from te head 1 detector (33.02 cm away from head 2). Therefore, their depths in water atthe two heads are 12.92 and 10.58 cm, respectively. In linear profile, MP-CFD generated result sketches the boundary of experimental data, and the log profile shows the difference between the scattered parts. From the comparison of the graphics in these two rows, it can also be seen that our result better matches the head 2 detection.
The accuracy of MP-CFD has also been investigated by simulating a SPECT acquisition of the NCAT phantom. We obtained noise-free projection images by HEGP collimators around the 64 × 64 × 64 NCAT phantom containing 100MBq of I-131. The accuracy of CFD has been verified in the previous validation part of this study by de Jong et al. and our previous work –. Therefore, we use CFD projection images to be the reference images. A total of 500 million photons were simulated for each CFD projection, while the MP-CFD required only a total of 500 million photon simulations for all 120 projection angles. Three base projection detectors are chosen (at the angle of α = 0°, α = 80°, and α = 200°). The NMSE value, as a function of detector angle, is introduced to evaluate MP-CFD using CFD resultant images to be reference images. It is plotted in Fig. 7. The maximum NMSE error is less than 0.03 for all measurements using different base projections. Fig. 8 depicts the comparison of resultant images of α = 0° base projection MP-CFD and CFD for different detectors when θ = 0°, θ = 120°, and θ = 240°. Several line profiles at the maximum error and other locations are analyzed, thereby verifying the accuracy of MP-CFD compared to CFD. The projection images of MP-CFD and CFD are well matched.
The same 64 × 64 × 64 NCAT phantom has also been used in speed tests of the MP-CFD method. A varying number of detectors (Np) have been simulated using MP-CFD, and the resultant projections are compared to CFD simulated projections with the same number of projections. For each simulation, total computation time for CFD and MP-CFD was recorded. As seen in Fig. 9, the simulation time of one CFD (MP-CFD) projection view is about 80 s, which includes 68 s photon propagation inside the phantom and 12 s of the attenuation and convolution calculation. With the increase ion the amount of the detector, the MP-CFD computation time increases because of the calculation of attenuation and convolution. The computation time is linear with the number of projections for both CFD and MP-CFD. The MP-CFD method is about six times faster than CFD.
We have previously assessed that the computation time of CFD is roughly 10 times faster than FD . This is due to the fact that the program uses an analytically determined blurring kernel rather than sampling the photon path direction from a probability density function. As a result, CFD converges to a noise-free projection image much faster than FD. The speed of MP-CFD is seen to be a further six times faster than CFD because the computationally demanding simulation of photon transport within the object only needs to be simulated once in MP-CFD and then forced to multiple detectors. Overall, the speed of MP-CFD is about 60 times faster than that of standard FD with similar image accuracy.
Contrary to the work of de Jong et al., photons are directed to the detectors at each photon interaction in our code, which results in the calculation of attenuation and convolution for all the photon transverse location. The attenuation coefficient is calculated with specific μ’s of each photon energy rather than the average , and the sub-projection map is no longer required to store the photon histories. Furthermore, the simulation is independent of the base projection because the combination of (4) and (6) cancels the attenuation coefficient of the base projection, which makes the photon history weights to other detectors only depend on the paths directed to them and the previous photon transmission weight but not the location of base projection.
In this paper, we have developed a practical method to accelerate Monte Carlo simulation of photon transport based on convolution-based forced detection (CFD). Following the simulated transport of photons in the object, the emitted photons are forced to be incident upon multiple detectors around the object simultaneously. This way, the computationally intensive task of modeling the photon transport through the object only needs to be done once for each photon. Furthermore, by analytically modeling the point spread function of the gamma camera and forcing the photon to interact within this entire region, additional speed improvements are made possible. With both multiple projection sampling and convolution-based forced detection, we have seen overall simulation times decrease by a factor of 60 compared to conventional forced detection (FD) while still maintaining similar accuracy.
This work was supported by the National Cancer Institute under Grant R01-CA104907-01A1. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the NCI.
Shaoying Liu, Department of Electrical and Computer Engineering, McMaster University, Hamilton, ON L8S 4K1 Canada.
Michael A. King, Division of Nuclear Medicine, University of Massachusetts, Medical School, Worcester, MA 01655 USA.
Aaron B. Brill, Department of Radiology, Vanderbilt, University, Nashville, TN 37232 USA.
Michael G. Stabin, Department of Radiology, Vanderbilt, University, Nashville, TN 37232 USA.
Troy H. Farncombe, Department of Nuclear Medicine, Hamilton Health Sciences, Hamilton, ON L8N 3Z5 Canada.