Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2931277

Formats

Article sections

Authors

Related links

IEEE Trans Nucl Sci. Author manuscript; available in PMC 2010 September 1.

Published in final edited form as:

IEEE Trans Nucl Sci. 2008 June; 55(3): 967–974.

doi: 10.1109/TNS.2008.924079PMCID: PMC2931277

NIHMSID: NIHMS230608

Shaoying Liu, Department of Electrical and Computer Engineering, McMaster University, Hamilton, ON L8S 4L8, Canada;

Shaoying Liu: ac.retsamcm@22suil

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) [1], [2]. 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 [3]. 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 [4].

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 [5]. Du, *et al.* used Monte Carlo N-Particle (MCNP) to simulate photon interaction types in the collimator [6]. 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 [7]. 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 [8]. 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 [3], [4], [9]. 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*(*O*_{1}*,*…, *O _{i},* …,

$${I}_{B}={I}_{0}exp(-{P}_{OB}\xb7\mu )$$

(1)

where, *I*_{0} 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.)

Representation of a point source above a parallel hole collimator. (a) Rays emitted from point source to the collimator. (b) Simplified illustration of the ray tracing method. (c) Collimator side view. (d) Collimator top view.

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
${A}^{\prime}({A}_{1}^{\prime},{A}_{2}^{\prime},\dots ,{A}_{n}^{\prime})$, 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 *P _{OB}* can be obtained by the calculation of the two dimensional (2D) path length

$${P}_{OB}=\frac{{P}_{{A}^{\prime}B}}{sin\theta}=\frac{{P}_{{O}^{\u2033}B}-{P}_{{O}^{\u2033}{A}^{\prime}}}{sin\theta}$$

(2)

where, *P _{O}*

If the location of B is denoted as (*L _{O}*

$$\theta =arctan\left(\frac{{L}_{{O}^{\u2033}B}}{z+l}\right)$$

(3)

where, the *L* represents the length of the segment denoted by the appropriate subscript, then *L _{o}*

$${L}_{{O}^{\u2033}{A}^{\prime}}=z\xb7tan(\theta )$$

(4)

and so the position of *A*′ can be written as (*L _{O}*

From the above discussion, it is easy to see that the intensity of photons at point *B* can be calculated by (1) using *P _{OB}* once the locations of

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*(*D*_{1}*,*…, *D _{i},* …,

Septal thickness calculation on the collimator back plane. (a) Rays on collimator back plane. (b) Schematic for checking whether point D is inside or outside a hexagonal hole.

After *P _{O}*

The 360° circle centered at point *O*″ is equally divided into *N _{r}* evenly spaced rays, as illustrated in Fig. 2(a). The location of points

Suppose the center of a hexagonal hole is given by location (*C _{x}, C_{y}*) with width and height denoted as

$$\{\begin{array}{l}\mid \stackrel{\sim}{y}\mid \phantom{\rule{0.16667em}{0ex}}<{\scriptstyle \frac{{L}_{y}}{2}}\hfill \\ \mid \stackrel{\sim}{y}-\sqrt{3}\stackrel{\sim}{x}\mid <{L}_{y}\hfill \\ \mid \stackrel{\sim}{y}+\sqrt{3}\stackrel{\sim}{x}\mid <{L}_{y}\hfill \end{array}$$

(5)

where, = *x* − *C _{x}*,

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:

$${S}_{y}={L}_{s}+{L}_{y}$$

(6)

and the center distance between 2 adjacent horizontal holes along a single row is:

$${S}_{x}=\sqrt{3}\xb7{S}_{y}$$

(7)

Equation (6) and (7) should be modified slightly for foil collimators. Assuming now that the current hole is located on the 0* ^{th}* row and 0

$$\{\begin{array}{l}{C}_{x}={\scriptstyle \frac{{S}_{x}}{2}}\xb7m\hfill \\ {C}_{y}={\scriptstyle \frac{{S}_{y}}{2}}\xb7n\hfill \end{array}$$

(8)

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 *L _{y}* <

The RT process starts with point *D* moving from *O*″ along one of the rays, with center hole (*C _{x}, C_{y}*) as (0,0) and initial cumulative septal thickness (

*D*_{i}_{−1}*H*and_{c},*D*_{i}*H*. This means point_{c}*D*is still within the same hole and thus_{i}*CST*=_{i}*CST*_{i}_{−1}. The current hole is still the center hole.*D*_{i}_{−1}*H*and_{c},*D*_{i}*S.*This means*D*is moving from air to septa._{i}*CST*=_{i}*CST*_{i}_{−1}+*d.*The current hole changes from the center hole to a side hole.*D*_{i}_{−1}*S,*the coordinates of*D*are required to be calculated with respect to each side hole using (5) to check whether it is inside one of the side holes. If_{i}*D*_{i}*S, CST*=_{i}*CST*_{i}_{−1}+*d,*and the current hole is still a side hole.*D*_{i}_{−1}*S,*but*D*_{i}*H*. This means_{s}*D*is moving from the septa to a collimator hole, and_{i}*CST*=_{i}*CST*_{i}_{−1}. The center hole is renewed as the one which holds the current*D*. The current hole changes back to the center hole, and the side holes are updated._{i}

This process is repeated until *D _{i}* goes to the end of all the rays being simulated with the

Once generated, the collimator response can be calculated using the table, *LUT _{PL},* by reprojecting the cumulative 2D septal thickness back to the 3D septal thickness based upon (2). By reading

$$\{\begin{array}{l}{V}_{i,j}={V}_{i,j}+{I}_{B}\hfill \\ {\mathit{CNT}}_{i,j}=CN{T}_{i,j}+1\hfill \end{array}$$

(9)

where, *i, j* denote the *i ^{th}* row,

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:

$${V}_{i,j}=\frac{{V}_{i,j}}{CN{T}_{i,j}}.$$

(10)

In order to generate the *LUT _{CR},* the point source is moved along the axis

The sampling intervals for the *d _{p}* and

Following the generation of *LUT _{CR},* it is then implemented into a CFD version of the SIMIND MC program (RT-MC) [1], [2]. In CFD, when a photon of energy

$$P={\mathit{PHW}}_{i}\times KN(\phi )\times (1-{P}_{pe}(E))exp(-\mu (E)\xb7{l}_{FP})$$

(11)

where, *i* is the photon interaction number, is the scattered angle, *l _{FP}* is the photon forced path, and exp(−

*LUT _{CP}* 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 (*L _{y}* = 0.4 cm,

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.

Comparison between RT-MC model and experimental result using HEGP collimator when point source located at 30 cm above the collimator. (a) and (b) compare the experimental measured data and RT-MC resultant image, which has been separated into geometric **...**

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 *LUT _{PL}*. Fig. 4(a) plots

The cumulative septal thicknesses from *O*″ and *O* to point *B* along a ray. The *x*-axis denotes the location of point *B* away from the center *O*″ with the unit of cm. (a) Cumulative septal thickness *P*_{O}_{″}_{B}. (b) The magnified view of the **...**

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 (*r*^{2}) has then been introduced to determine how well the RT models fit the experimental data.

RT method sensitivity evaluation for PSF with point source in air (left) and in water (right). It presented the sensitivity of the detection based upon the detected photon counts per second for 1MBq activity.

Table I shows the *r*^{2} evaluation of HEGP collimator for sources in both air and water, together with the same evaluations of medium energy general purpose collimator (MEGP collimator: *L _{y}* = 0.3 cm,

We can see that all the *r*^{2} 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 [10]. 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 [5].

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).

cv vs time. The solid line is the computation time of rt, and the dashed line is computation time of non-RT mc. The computation time of rt-mc is much less that non-RT mc, and the stable noise level of it is lower than non-RT mc.

In the implementation of the RT method, it is possible to alter the number of the rays (*N _{r}*), ray step size (

An optimization has been performed to evaluate the minimum *N _{r}* necessary, and the optimal

$$\mathit{NMSE}=\frac{\sum _{i=1}^{n}{({O}_{i}-{F}_{i})}^{2}}{\sum _{i=1}^{n}{{F}_{i}}^{2}}$$

(12)

where, *O _{i}* and

Fig. 9(a) shows the average value of the NMSE for the same *N _{r}* when the point source moves from 10 cm to 60 cm in 10 cm intervals, and

NMSE comparison of the model for the 3 collimators with the reference projection images over various *N*_{r} and *d.* (a) shows the effect of *N*_{r} is not very obvious, 1,000 rays is enough for simulation. (b) shows the effect of *d* is significant. The optimal value **...**

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 *LUT _{CR}* 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,

For the determination of *d _{p},* we have simulated the PSF’s for two point sources with a specified

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.

1. de Jong HWAM, Slijpen ETP, Beekman FJ. Acceleration of Monte Carlo SPECT simulation using convolution-based forced detection. IEEE Trans Nucl Sci. 2001 Feb;48(1):58–64.

2. Farncombe TH, Liu S, King MA, Brill AB, Stabin MG. Accelerated SPECT Monte Carlo simulation using convolution-based forced detection. Proc Soc Nucl Med Annu Conf. 2006

3. Beck RN, Redtung LD. Collimator design using ray-tracing techniques. IEEE Trans Nucl Sci. 1985 Feb;NS-32(1):257–272.

4. Muehllehner G, Luig H. Septal penetration in scintillation camera collimators. Phys Med Biol. 1973;18(6):855–862. [PubMed]

5. Ljungberg M, Larsson A, Johansson L. A new collimator simulation in SIMIND based on the delta-scattering technique. IEEE Trans Nucl Sci. 2005 Oct;52(5):1370–1375.

6. Du Y, Frey EC, Wang WT, Tocharoenchai C, Baird WH, Tsui BMW. Combination of MCNP and simSET for Monte Carlo simulation of SPECT with medium- and high-energy photons. IEEE Trans Nucl Sci. 2002 Jun;49(3):668–674.

7. Song X, Segars WP, Du Y, Tsui BMW, Frey EC. Fast modelling of the collimator-detector response in Monte Carlo simulation of SPECT imaging using the angular response function. Phys Med Biol. 2005;50:1791–1804. [PubMed]

8. Staelens S, de Wit T, Beekman F. Fast hybrid SPECT simulation including efficient septal penetration modelling (SP-PSF) Phys Med Biol. 2007;52:3027–3043. [PubMed]

9. Lewis DP, Tsui BMW, Tocharoenchai C, Frey EC. Characterization of medium and high energy collimators using raytracing and Monte Carlo methods. Proc. 1998 IEEE Nucl. Sci. Symp. Med. Imaging Conf.; 1998. pp. 2026–2030.

10. Gunter DL. Collimator design for nuclear medicine. In: Wernick MN, Aarsvold JN, editors. Emission Tomography-The Fundamentals of PET and SPECT. New York: Academic; 2004. pp. 153–168.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |