If you have trouble viewing this, try the pdf of this post. You can download the code used to produce the figures in this post.

In this post, I continue the discussion of my paper “*Near optimal energy selective x-ray imaging system performance with simple detectors*[1]”, which is available for free download here. I will describe a Monte Carlo simulation of the SNR with energy information discussed in my last post. The simulation traces the paths of individual photons. This is, perhaps, more fundamental than my previous simulations that relied on statistical models for the detector data. I develop models for the random path lengths and use them to simulate the imaging task for SNR described in my last post. The model is validated by comparing the energy spectrum of photons transmitted through an object to the theoretical formula. I also provide an estimate of the errors in the Monte Carlo results.

A good reference for this topic is the book *A Monte Carlo Primer* by S. A. Dupree and S. K. Fraley[2], although the book is very expensive. I cover all the necessary topics here so the book is not necessary to understand this post but it covers many topics that I do not discuss.

As discussed in my previous post, we can describe the probability of a photon interacting with matter by a cross section *σ*(*E*), which is defined by

In this equation, a thin sample of the material is placed between collimators so any interaction removes it from the beam. If *n* photons of energy *E* are incident on a material of thickness *δ**x* with *N* atoms per unit volume., then, on average, *n* − *δ**n* photons pass through the material without interacting. Integrating Eq. 1↑, the number that have not interacted after a thickness x will be

The linear attenuation coefficient is

which has the units inverse length *m*^{ − 1}.

Eq. 2↑ is the same as the probability distribution for the inter-pulse time of a Poisson process discussed here, where the “rate” is *μ*(*E*). Notice that the path length distribution depends on the photon energy. Fig. 1↓ shows the geometry and the imaging task. We measure the number of photons **n**_{b} transmitted through the background region and through a region with background and possibly a feature **n**_{f} and use the data to decide whether the feature is present.

In the simulation, photons **n**_{0} are generated with a specified energy distribution. From the energy *E* of an individual photon, the attenuation coefficient is calculated and used to generate a random path length

where *rand* is a uniform (0, 1) random number. The model I use is that any interaction with the material stops the photon so the only ones that get through the background region are those for which **t** > *t*_{b}. I also assume that the detector measures all photons that get through the object. These are obviously an oversimplification but they capture the essential features of the problem. In Fig. 1↓, photons 1 and 3 pass through the material but number 2 does not.

The feature region is more complex since it has two materials with different *μ*(*E*). There are two ways to approach this problem. One is to express the random path lengths in terms of mean free path (MFP). Then the random paths are simply **t** = − log(*rand*) and we can convert the MFP to the actual path length with a function. However the conversion depends on the photon energy so I decided to use an alternate approach. This is a two step process. The first step is the same as for the background-only region. We generate photons, compute their attenuation coefficient and use this to compute a random path length. If this is greater than the thickness of the background material, then we “re-start” the photon at the boundary between the two materials regardless of the actual path length. We compute a new path length using the attenuation coefficient of the feature material and if it is greater than *t*_{f} the photon is counted. This approach is valid because a Poisson process is “memory less.” The probability of interacting is constant regardless of the position along the path. Once a photon passes the intersection we can compute a new random path starting at the intersection and not affect the results.

Some example photon paths are shown in Fig. 1↓. Photons 4 and 5 are not counted. Although the initial path for photon 4 was greater than *t*_{b}, the random path length from the intersection is less than *t*_{f} so the photon does not pass through. The initial path length of photon 5 was less than *t*_{b} − *t*_{f}, so it is not counted. Photon 6 is counted because its initial path length is greater than *t*_{b} − *t*_{f} and the second path length is greater than than *t*_{f} .

The simulation uses two functions. The first one *MakePulses* is shown below. First it generates random inter-pulse times, *dtpulse* and their cumulative times, *ttotal*. The number of pulses used is the minimum so the cumulative time is greater than the specified integration time *t*_*integrate*. The function then generates random energies for these photons using the inverse cumulative transform method discussed in this post.

Based on the photon energies computed by *MakePulses*, the function *MakePathLengths* computes a set of random path lengths. The attenuation coefficient for each energy is computed using linear interpolation by the *interp*1function. These attenuation coefficients are then used to compute the random path lengths using Eq. 4↑.

The software was tested by computing the spectrum of transmitted photons. The relevant code is shown below

The theoretical formula for the spectrum of the transmitted photons is
*n*_{b}(*E*) = *n*_{0}(*E*)*e*^{ − μbtb}

Fig. 2↓ shows the histograms of the incident and the transmitted photon spectra and the theoretical formulas. The spectra are scaled so they have the same maximum as the histograms. although the transmitted spectra is noisy since there are many fewer photons than the incident spectrum, the Monte Carlo data match the theoretical formulas well.

A large number of trials were used to compare the SNR with different types of detectors. The SNR was computed as the difference of the mean values in the feature and background region squared divided by the variance of the background region data. The photons transmitted through the feature region were computed using the two step algorithm described above. The relevant code is shown below:

The counts weighted with the optimal function described in in my last post were computed as shown below. Since the photons are processed individually, this is just the sum of the weights for the photon energies.

% each photon occurrence is weighted once Nwgtft = sum(interp1(egys,wgts,egys_transmitted));

Fig. 3↓ is a plot of the SNR with each of the detector types. The Monte Carlo results validate the theoretical result
*SNR*^{2}_{optimal_wgt} ≥ *SNR*^{2}_{photon counting} ≥ *SNR*^{2}_{energy integrating}

Fig. 3↑ gives an idea of the gain in SNR with complete energy information. The optimal SNR is about twice the energy integrating detector, which is the present state-of-the-art. Since the *SNR*^{2} is inversely proportional to the dose, the optimal system can provide the same SNR with less than half the dose. That is a significant advantage and well worth pursuing.

The next step is to develop theory that allows us to compare the SNR with limited energy-resolution detectors to the optimal value and to each other. This will lead me to a discussion of statistical detection theory.

Last edited Nov 15, 2011

© 2011 by Aprend Technology and Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.

[1] : “Near optimal energy selective x-ray imaging system performance with simple detectors”, *Med. Phys.*, pp. 822—841, 2010.

[2] : *A Monte Carlo Primer: A Practical Approach to Radiation Transport*. Springer, 2001.