|Home | About | Journals | Submit | Contact Us | Français|
In SPECT imaging, photon transport effects such as scatter, attenuation and septal penetration can negatively affect the quality of the reconstructed image and the accuracy of quantitation estimation. As such, it is useful to model these effects as carefully as possible during the image reconstruction process. Many of these effects can be included in Monte Carlo (MC) based image reconstruction using convolution-based forced detection (CFD). With CFD Monte Carlo (CFD-MC), often only the geometric response of the collimator is modeled, thereby making the assumption that the collimator materials are thick enough to completely absorb photons. However, in order to retain high collimator sensitivity and high spatial resolution, it is required that the septa be as thin as possible, thus resulting in a significant amount of septal penetration for high energy radionuclides. A method for modeling the effects of both collimator septal penetration and geometric response using ray tracing (RT) techniques has been performed and included into a CFD-MC program. Two look-up tables are pre-calculated based on the specific collimator parameters and radionuclides, and subsequently incorporated into the SIMIND MC program. One table consists of the cumulative septal thickness between any point on the collimator and the center location of the collimator. The other table presents the resultant collimator response for a point source at different distances from the collimator and for various energies. A series of RT simulations have been compared to experimental data for different radionuclides and collimators. Results of the RT technique matches experimental data of collimator response very well, producing correlation coefficients higher than 0.995. Reasonable values of the parameters in the lookup table and computation speed are discussed in order to achieve high accuracy while using minimal storage space for the look-up tables. In order to achieve noise-free projection images from MC, it is seen that the inclusion of the RT implementation for septal penetration increases the speed of the simulation by a factor of about 7,500 compared to the conventional SIMIND MC program.
Monte Carlo (MC) is a well-known computational tool used to simulate the photon transport in single photon emission computed tomography (SPECT). While accurate, the main drawback of MC is its enormous computational burden which prevents it from being routinely implemented into clinical practice. In recent years, several variance reduction techniques (VRT) have been developed to improve its detection efficiency as well as the computation speed. One of the latest implementations used to significantly accelerated MC modeling is convolution-based forced detection (CFD) , . With CFD, photons are tracked through the object under study, and then forced to travel in the direction perpendicular to detector surface. The probability of photon detection is convolved with a Gaussian blur function representative of the collimator response function. This mathematical model can only estimate the geometric response based on the assumption that photons traveling in the collimator material will be completely attenuated. However, when the photon energy exceeds 200 keV, there is high probability that photons are able to penetrate the collimator . Furthermore, some photons may interact with the collimator materials before they penetrate the septa. Therefore, a complete characterization of the collimator response involves geometric response, septal penetration and scatter. Although different collimators are used to minimize septal penetration, thin septa are preferred when high detection efficiency is required. Therefore, the effect of septal penetration is typically significant when high or medium energy radionuclides are used. Overall, the amount of septal penetration is determined by the path length of rays through the septa which depends on the distance between the collimator and radioactive source, hole length, hole pattern and septal thickness .
Several studies have been done to model collimator response. Ljungberg, et al. have developed a Delta-Scattering technique by sampling a series of random numbers to determine the photon interaction types in the collimator when the traverse of photons are tracked . Du, et al. used Monte Carlo N-Particle (MCNP) to simulate photon interaction types in the collimator . Song, et al. presented a method using MC code to generate the collimator response models based upon the angular response functions (ARFs), which depends on the photon traveling directions and photon energies . Recently, Staelens, et al. extended convolution-based forced detection method (CFD) by convolving the photon transmission probabilities with MC pre-simulated collimator responses when ray tracing method is performed onto the object . However, the MC simulation of a noise-free collimator response image is extremely time consuming. In this work, we will introduce a method using ray tracing (RT) to analytically generate the collimator response model including geometric response and septal penetration. The collimator scatter compartment is not going to be discussed here because of the difficulty in analytically simulating photon scatters. RT has been previously used to evaluate the collimator characters by modeling the collimator response , , . However, it has not been applied to compensate for collimator response. Eventually, the RT resultant models will be applied into CFD MC code for further acceleration.
The fundamental idea of the RT method is shown in Fig. 1(a). Consider a parallel hole collimator and a point source located at the distance, z, above the collimator surface. The gray regions denote the collimator septa and the white represent the collimator hole. A ray passing through the gray region thus contributes to the septal penetration, otherwise, it contributes to the geometric response. Photons are simulated as emitting isotropically from the point source O(O1,…, Oi, …, On), thus striking onto the collimator back plane. Suppose one of the emitted rays interacts with the collimator at point A(A1, …, Ai, …, An) on the front plane and B(B1, …, Bi,, …, Bn) on the back plane, where the subscript i corresponds to ith ray. The direction of each ray can be specified by the polar angle γ(γ1, …, γi, …, γn) and azimuth angle θ( θ1, …, θi, …, θn). The intensity of the photons striking point B depends on the path length of OB in the septal region, as denoted by POB. It is easy to see that POB = 0 for geometric response (i.e., no septal penetration). The equation of intensity can be written as:
where, I0 is the initial photon interaction probability at a certain photon scatter location, and μ is the linear attenuation coefficient of the collimator material at a specific photon energy (e.g., 3.12 cm−1 for lead at 364.5 keV.)
Fig. 1(b) and (c) simplify Fig. 1(a) and specify two interesting ray characteristics which stress the fundamental ideas of the RT method as presented in this paper. The first point is that the projections of rays from the source along the axis OO″ to a certain point B on the collimator are the same no matter where the point source O is located. Thus, projection line of ray OB is always O″B when point source O is moved along the axis OO″. In addition, the location of point , which is the projection of point A, is moved along O″B with respect to the movement of O, as shown in Fig. 1(b). The second point is that the three dimensional (3D) septal path length POB can be obtained by the calculation of the two dimensional (2D) path length PA′B, which denotes the path length of A′B in the gray region. We then have:
where, PO″B and Po″A′ represent the path length of O″B and O″ A′ in the gray region, respectively. It is noted that (2) can be only used for the case of a parallel hole collimator.
If the location of B is denoted as (LO″B,γ) (LO″B is the length of segment O″B, and γ is the angle of O″B away from the horizontal direction.), and defined by polar coordinates as shown in Fig. 1(d), the value of θ and location of point A′ in Fig. 1 can be derived from the collimator parameters and location of the point source O. Suppose the point source is placed above collimator at a distance z, the collimator length is l, and the angle of ray O″ A′ to the horizontal direction is γ, the same as ray O″B. θ can be derived from the trigonometric functions:
where, the L represents the length of the segment denoted by the appropriate subscript, then Lo″A′ is:
and so the position of A′ can be written as (LO″A′ γ) on polar plane.
From the above discussion, it is easy to see that the intensity of photons at point B can be calculated by (1) using POB once the locations of B and O are known. For different O positions, the corresponding values of PO″B are the same while PO″A″ may not. However, PO″A′ may be the same for different B. The calculation of PO″B and PO″A′ are significantly time consuming, therefore, it would be helpful to produce a lookup table LUTPL to store the cumulative septal thickness PO″B and PO″A. In this study, hexagonal hole parallel collimators are used, although this method is applicable for all types of parallel collimators.
Several rays are traced in the collimator back plane from one center point O″ in different directions, with each ray modeled by a large number of equally spaced segments denoted D(D1,…, Di, …, Dn), as seen from Fig. 2(a). A region including a collimator hole and the surrounding septa is equally divided into n sub-regions, and the center location of each sub-region is set as the center point O″. n is a positive integer number which should at least be greater than the ratio of the hole-septa area to a projection image pixel area. In the following discussion, O″ is chosen as the center of a collimator hole as an example.
After PO″D, the septal path lengths of segment O″D are calculated and stored in LUTPL together with their polar coordinates, a second lookup table consists of the resultant collimator response functions (denoted as LUTCR) based upon LUTPL and is constructed by the projection images of a point source located at different distances from the collimator at a variety of photon energies. LUTPL is incorporated into SIMIND MC program in order to accelerate it. The Section II-B discusses the creation of LUTPL and LUTCR and their subsequent implementation in CFD-MC.
The 360° circle centered at point O″ is equally divided into Nr evenly spaced rays, as illustrated in Fig. 2(a). The location of points D on the rays are computed to see whether they are within the septa or inside the collimator hole, which can be seen in Fig. 2(b).
Suppose the center of a hexagonal hole is given by location (Cx, Cy) with width and height denoted as Lx and Ly, respectively, and the distance between two adjacent points Di−1 and Di is denoted as d. d should be satisfied with the conditions that d Ls and d Ly to minimize the error. A study of the amount of rays Nr and sampling interval d is performed below in the Discussion. A given point (x, y) within the hexagonal hole should satisfy:
where, = x − Cx, ỹ = y − Cy. These three inequality expressions correspond to the functions of the 6 edges of the hexagonal hole. Note this equation varies with different hole geometries.
When point D is moved along a ray, it is difficult to determine in which hexagonal hole the point D is contained. With so many holes in one collimator, it is exceedingly time consuming to determine all the hexagons by (5). Rather, suppose the coordinates of point O″ is (0,0), then the center distance between two adjacent vertical holes along a single column is:
and the center distance between 2 adjacent horizontal holes along a single row is:
Equation (6) and (7) should be modified slightly for foil collimators. Assuming now that the current hole is located on the 0th row and 0th column, (i.e., location (0,0)), then the coordinates for the hole centered on the nth row and mth column are:
where n and m can be any integer number.
Observing the hole pattern and septal thickness, when a ray exits out of a hole and proceeds to the next hole, the second hole is always adjacent to the first. For a practical (low/medium/high energy, general purpose/high resolution) collimator, it is usually adequate to consider only the adjacent 18 holes, as depicted by the shaded holes in Fig. 2(a). The number of adjacent holes should be increased only when 2 Ly < Ls, which makes the next point D always traveling in the septa until it gets into a new hole other than the shaded hole. Initially, we will call the current hole that D resides, the center hole, and its 18 adjacent holes as side holes.
The RT process starts with point D moving from O″ along one of the rays, with center hole (Cx, Cy) as (0,0) and initial cumulative septal thickness (CST0) as 0 cm, with the subscript 0 denoting the starting D. The previous point D is recorded as Di−1, and the current one being Di. When D is inside the center hole (one of the side holes), it is denoted as D Hc (D Hs), with D within the septa being denoted as D S. During the movement along ray O″D, four possibilities may occur:
This process is repeated until Di goes to the end of all the rays being simulated with the CST of each position on the ray being stored in the lookup table LUTPL.
Once generated, the collimator response can be calculated using the table, LUTPL, by reprojecting the cumulative 2D septal thickness back to the 3D septal thickness based upon (2). By reading PO″ B in Fig. 1 from LUTPL, the polar coordinates (LO″ B,γ) of B and the cumulative septal thickness PO″ B are known. The coordinates (LO″ A′, γ) of A′ can be derived from (3) and (4), and the subsequent value of PO″ A′ can be determined by searching LUTPL Furthermore, POB can be determined based on (2). Therefore, the intensity of the photon flux striking detector location B is determined by (1). As the detector surface is typically divided into a number of detector bins, this intensity gets added to the value of the detector bin corresponding to this location and the number of photons striking at this bin is subsequently incremented.
where, i, j denote the ith row, jth collimator bin. Vi,j is the value of this bin, IB is the intensity of certain photon flux on point B, as calculated from (1). CNTi,j is the total amount of the photons striking on this bin.
This procedure is repeated for all the locations on a given ray and for all possible emission angles over 360°. Eventually, the value of each bin can be given by:
In order to generate the LUTCR, the point source is moved along the axis OO″ at a step of dp and the value of μ at a step of dμ in (1) is altered according to different photon energies. The sum of each resultant projection image is normalized to 1 before the image is stored in LUTCR for each point source location and different μ. The matrix size of RT resultant images should be twice of the projection images to ensure the corners of the projection images get effective data when the photons are traveling above the opposite corners on the diagonals.
The sampling intervals for the dp and dμ are analyzed in the Discussion.
Following the generation of LUTCR, it is then implemented into a CFD version of the SIMIND MC program (RT-MC) , . In CFD, when a photon of energy E, undergoes scatter within the object, it is directed to the path perpendicular to the detector. The photon probability P is determined by (11), which is the product of photon history weight (PHW) beginning with 1 at the initial emission location, the probability of Compton scattering, as calculated by the Klein-Nishina (KN) function, and the probability (1 − Ppe(E)) that the photon is not absorbed by photoelectric interaction.
where, i is the photon interaction number, is the scattered angle, lFP is the photon forced path, and exp(−μ(E)·lFP) is the attenuation along the forced path.
LUTCP is then consulted to determine the projection images that most likely matches the photon energy and its scatter location. If an exact match is not found, then the projection image is obtained via linear interpolation.
In order to evaluate RT-MC generated models, a GE Millennium VG camera is simulated and compared with experimental data for a block water phantom. A high energy collimator (Ly = 0.4 cm, Ls = 0.18 cm) and I-131 point sources are used with 10,000 rays simulated over 360°, and with d is set to 0.004 cm. A 10% energy window centered at 364.5 keV is used to create the projection images.
Fig. 3 shows the results of RT-MC along with experimental data for I-131 when the point source is placed at 30 cm above the collimator in the air. The RT-MC generated image matches the experimental data very well. Additionally, the RT image for the HEGP collimator has been rotated counterclockwise by 15° to match experimental results. The center dark circle region corresponds to the geometric response with the long tails of each image corresponding to septal penetration, which have been individually illustrated in Fig. 3(c) and (d). Septal penetration appears as very well-known 6-sided star pattern with a second, less obvious star pattern situated 30° away from the first one. This can be explained in Fig. 1(d), we can see that in the vertical direction (black lines), the rays at every 60° travel through the thinnest septal thickness. The second, less intense, star pattern appears offset 30° (gray lines) from the first pattern due to another thin septal thickness compared to its surrounding direction. Fig. 3(e) and (f) respectively give the linear-plot and log-plot of one tail profile to show the difference between RT-MC and experimental data. RT-MC simulation shows more greatly fluctuated ripples is because our model is lack of the photon interactions such as collimator scatter, X-ray emission, and else interactions after crystal.
There is an additional structure seen in the PSF images as regions of higher and lower intensity along the septal penetration tails in Fig. 3. This can be explained from the CST’s stored in the LUTPL. Fig. 4(a) plots PO″B in the horizontal direction for a single ray. Fig. 4(b) shows a magnified view of the selected region in (a), and Fig. 4(c) and (d) plot the septal thickness PA′B and PAB for differing B when the point source is located 30cm above the collimator, CST appears as a staircase pattern, with the flat lines corresponding with point B in the collimator hole and the increasing lines associated with the points in the collimator septa. In Fig. 4(c) and (d), we can see some ripples in this pattern whereby several ripples combine to a larger packet. These ripples and bumps come from the numerator term in (2). The increasing portions or flat portions in Fig. 4(a) minus another increasing portion results in periodically appearing ripples. Between every 2 bumps, PO″B and PO″A moves along two flat lines and the minus term results in the slowly increasing line, which denotes the increasing azimuth angle canceling with the increase in the septal thickness. The sin term in the denominator of (2) flattens the septal thickness PAB, but enlarges the first bump more than others because of the small photon emission angle.
Full width at half maximum (FWHM) has been introduced to evaluate the RT-MC geometric compartment. Because the result of RT does not include the intrinsic resolution of the detector, it is necessary to further convolve the result of RT with the camera intrinsic resolution. This has been found to consist of a σ = 0.2 cm Gaussian function for HEGP collimator using 2 cm thick crystal. The experimental and RT-MC generated images are convolved with an additional σ = 0.2 cm Gaussian function, to remove the effect of the collimator hole pattern as seen in Fig. 3. Fig. 5 shows the FWHM evaluation of RT and experimental PSF’s when the point source is located in the air and water, respectively. The correlation coefficient (r2) has then been introduced to determine how well the RT models fit the experimental data.
Table I shows the r2 evaluation of HEGP collimator for sources in both air and water, together with the same evaluations of medium energy general purpose collimator (MEGP collimator: Ly = 0.3 cm, Ls = 0.1 cm) with In-111 and low energy high resolution collimator (LEHR collimator:Ly = 0.15 cm, Ls = 0.02 cm) with Tc-99m when the source is placed in the air.
We can see that all the r2 of the PSF’s of these three collimators are greater than 0.995, showing that the RT models are consistent with experimental data. Fig. 5(a) and (b) shows that the simulated detector sensitivity (CPS/MBq) of RT-MC agrees well with the measured data, for both air or water.
Full width at tenth maximum (FWTM) is taken as the measurement criteria to assess the accuracy of the septal penetration model . Excellent agreement between the RT-MC models and the experimental PSF’s for these three collimators is also seen in the Table I.
RT models do not include the photon interactions in the collimator and gamma camera. However, photon scatter within the collimator and crystal, X-ray emission, and backscatter from the materials behind the crystal can be very important when medium or high energy photons are simulated. Fig. 6 shows the energy spectra of I-131 using HEGP collimator. Photon interactions within the collimator and gamma camera affect the photons with energies lower than 200 keV much more than those with higher energies. In practice, the 10% energy window centered at 364.5 keV is used, which makes the photon interactions not a significant problem. For more accurate simulation, the new SIMIND MC code developed by Ljungberg, et al., is recommended to be combined with the RT models .
In order to evaluate the speed of the RT method when implemented into the SIMIND MC program, we have simulated a 50 × 50 × 50 I-131 block source contained within a 128× 128 × 128 voxelized air space using a HEGP collimator. Different numbers of photons have been simulated and the computation time recorded. As the projection image is expected to be relatively uniform, the coefficient of variation (CV) of the center 50 × 50 dark region in the resulting images is calculated to measure the degree of uniformity.
Fig. 7 shows a plot of CV vs computation time. As seen, the RT method requires approximately 100 seconds in order to achieve a stable level of CV and a projection image with CV = 0.006. Forced detection MC is performed as a non-RT Monte Carlo method (non-RT MC) to compare with RT-MC. It takes the non-RT SIMIND program 1494 seconds to achieve an image with CV = 0.03. In order to achieve a noise free image with CV = 0.006, it takes the non-RT MC longer than 750,000 s (9 days) to run, indicating a speed improvement for the RT-MC of about 7,500 times. Fig. 8(a) and (b) shows the projection images of non-RT MC and RT-MC with CV = 0.006. It is easy to see that the result of RT-MC is relatively noise free but even after 9 days simulation, the non-RT MC is not. The center regions both depict the hexagonal pattern, which is due to the geometry of the collimator. The resultant non-RT image looks sharper than RT possibly because the maximum emission angle of septal peneration rays is limited in SIMIND MC for difference source-detector distances. Only one single maximum emission angle can be used for the cubic source which include several different distances. The limitation will result in insufficient septal penetration photons and an overestimate of the geometric photons, which can also been seen in the linear and log center profiles shown in Fig. 8(c) and (d).
In the implementation of the RT method, it is possible to alter the number of the rays (Nr), ray step size (d), the distance between two adjacent source/detector planes (dp) and the linear attenuation coefficient spacing (dμ) used in the LUTPL corresponding to different photon energies and LUTCR lookup tables. By increasing Nr, and by reducing d, dp, and dμ, the resultant septal penetration model should yield higher accuracy, but at the expense of increased computation time and storage requirements for LUTPL and LUTCR. As such, a study of these parameters has been undertaken to determine the optimal setting in order to maintain accuracy, yet reduce computational and storage requirements. The following section discusses the appropriate parameters for HEGP, MEGP, and LEHR collimators.
An optimization has been performed to evaluate the minimum Nr necessary, and the optimal d in order to minimize computer requirements while maintaining accuracy. The CST and Ly of different collimators vary, which makes the optimal simulated parameters of each collimator different. The evaluation depends on the hole size and cumulative septal thickness. For reference, we have used projection images generated using a lookup table calculated with Nr = 10,000 and d = 0.002 cm, which have shown high accuracy compared with experimental data. PSF’s modeled using reduced numbers of rays, Nr, and larger step sizes, d, have been compared to these reference images using the normalized mean square error (NMSE):
where, Oi and Fi are the ith pixel values of the estimated projection images and reference images, respectively.
Fig. 9(a) shows the average value of the NMSE for the same Nr when the point source moves from 10 cm to 60 cm in 10 cm intervals, and d is altered from 0.002 cm to 0.15 cm. Fig. 9(b) shows the average value of the NMSE for the same d when the point source moves from 10 cm to 60 cm at each step of 10 cm, and Nr is altered from 10,000 to 500. Fig. 9(c) and (d) are the log value of Fig. 9(a) and (b). We can see that d affects the result more than Nr. In fact, it can be seen that the LEHR collimator requires more precise parameters due to its thinner septal thickness and smaller hole size. Table II shows the optimal value of Nr and d when the corresponding NMSE value gets to a stable level. For all collimators, it is seen that 1,000 rays are enough for RT simulation when the step size is set to half of the collimator septal thickness.
As the location of an emitted photon is variable and the resultant collimator attenuation coefficient varies widely as a function of energy (ranging from 1.167 cm−1 at 700 keV to 62.13 cm−1 at 100 keV), it is impossible to store all possible source location and energy combinations in the lookup table. Rather, we have implemented the LUTCR for set distances from the detector surface, and for set photon energies. Thus, for a given source location and photon energy, we determine the correct collimator response model via linear interpolation between the stored values. Thus, we must determine the optimal spacing for source-detector distance, dp and linear attenuation coefficient, dμ.
For the determination of dp, we have simulated the PSF’s for two point sources with a specified dp. A linear interpolation is performed between these PSF’s and compared to the actual PSF at that location using the NMSE. The evaluation for dμ is performed in an analogous manor. HEGP collimator are used in the following evaluation. As seen in Fig. 10, when the separation between calculated values for dμ increases, the NMSE increases as the linear interpolation scheme is not accurate. It can generally be seen that reasonable accuracy exists in the simulated PSF even for dμ = 0.5 cm−1. As expected, it is seen that by reducing the plane separation, dp, more accurate PSF’s are obtained. This trend is seen to be essentially linear, thereby suggesting that the optimal plane spacing may be equal to the voxel size as simulated in the subsequent MC simulation.
In this work, we have introduced an effective ray tracing code into the SIMIND Monte Carlo program to account for septal penetration within the collimator. In order to efficiently implement this method, a lookup table storing the cumulative septal path lengths for a point source must first be generated. Once evaluated and stored, specific collimator response functions can then be calculated based on the path length lookup table. The collimator response functions are then stored in a subsequent lookup table for different source–detector distances and various photon energies. In order to accelerate Monte Carlo modeling, the collimator response lookup table is referenced for each photon that is emitted from the object simulated. The point spread function referenced from the lookup table is then convolved with the detection probability on the camera surface, thereby resulting in convolution-based forced detection incorporating septal penetration modeling.
An optimization of various parameters stored in the lookup tables has been performed for different types of collimators. A validation has been performed assessing the correlation coefficient for the point spread functions, and has shown that the ray tracing method is very accurate. When incorporated into the SI-MIND Monte Carlo program, it has been shown that the speed of the ray tracing method has increased the simulation speed by a factor of 7,500 compared to the conventional Monte Carlo program. Subsequent work will focus on using this RT-MC for image reconstruction in SPECT imaging.
Shaoying Liu, Department of Electrical and Computer Engineering, McMaster University, Hamilton, ON L8S 4L8, Canada.
Michael A. King, Department of Radiology, University of Massachusetts Medical School, Worcester, MA 01655 USA.
Aaron B. Brill, Department of Radiology, Vanderbilt University, Nashville, TN 37240 USA.
Michael G. Stabin, Department of Radiology, Vanderbilt University, Nashville, TN 37240 USA.
Troy H. Farncombe, Department of Nuclear Medicine, Hamilton Health Sciences, Hamilton, ON L8R 1M8, Canada.