According to Equation 6.2, the final expected number of hits is given by a long weighted sum of the expected number of events between the neighboring crystals, i.e. the LOR value obtained with the geometric model. Note that this is similar to image filtering, but now the space is not 2D but 4D (see Section 6.2). The long sum is evaluated by MC estimation taking Ndet random samples of detector pairs (i(1),j(1)),(i(2),j(2)), . . . ,(i(Ndet),j(Ndet)):
˜
ydetL(d1,d2)≈µd1µd2C(˜ygeom,d1,d2) =∑
i
∑
j
˜
yL(i,j)geom·pi→d1 ·µd1 ·pj→d2 ·µd2 ≈
1 Ndet ·
N∑det
s=1
˜
yL(i(s),j(s))geom ·pi(s)→d1·µd1 ·pj(s)→d2 ·µd2
ps
wherepsis the probability of samples. A sample is associated with a pair of offset vectorsd1−i from d2−j. According to importance sampling [Chr03], ps is made proportional to the crystal transport probability:
ps= pi(s)→d1 ·pj(s)→d2
∑
i
∑
jpi→d1 ·pj→d2 = pi(s)→d1 ·pj(s)→d2 ν1(⃗ω)·ν2(⃗ω) . Thus, the final estimator is:
˜
yL(ddet1,d2)≈ ν1(⃗ω)·ν2(⃗ω)·µd1 ·µd2
Ndet ·
N∑det
s=1
˜
ygeomL(i(s),j(s)). (6.5) This method runs a geometric first pass, which is the same algorithm as developed to execute the forward projection of the geometric reconstruction. This pass results in LOR values ˜ygeomL . Then, the 4D LOR map is filtered. We visit again each LOR, find neighbors of its two crystals according to a prepared random map, and add up the values stored in the LOR selected by the two sampled neighbors.
CHAPTER 6. DETECTOR MODEL 65 6.3.1 Detector modeling in back projection
In back projection the system matrix is simplified and we ignore blurring effects like positron range and scattering. Using piece-wise constant approximation of the transport function we get:
AL(d1,d2),V ≈∑
i
∑
j
DL(i,j),V ·pi→d1(⃗ωd1→d2)·µd1·pj→d2(⃗ωd1→d2)·µd2 whereDLV is the system matrix simulating geometric effects and attenuation.
In the back projection of ML-EM reconstruction, we have to evaluate the numerator and the denominator of the scaling factor for each voxel. The numerator is
∑
L
ALV ·yL
˜
yL =∑
d1
∑
d2
AL(d1,d2),V ·yL(d1,d2)
˜
yL(d1,d2). Let us substitute the factorization of the system matrix into this expression:
∑
L
ALV ·yL
˜
yL ≈∑
d1
∑
d2
∑
i
∑
j
DL(i,j),V ·pi→d1(⃗ωd1→d2)·µd1·pj→d2(⃗ωd1→d2)·µd2 ·yL(d1,d2)
˜
yL(d1,d2) =
∑
i
∑
j
DL(i,j),V ·∑
d1
∑
d2
pi→d1(⃗ωd1→d2)·µd1·pj→d2(⃗ωd1→d2)·µd2 ·yL(d1,d2)
˜
yL(d1,d2).
As the crystal transport probability depends just on the distance, we can reverse the direc-tion:
pi→d1(⃗ω) =pd1→i(−⃗ω).
Note that the nominator can also be expressed as a geometric back projection from a term obtained with convolution:
∑
L
ALV ·yL
˜
yL ≈∑
i
∑
j
DL(i,j),V · C
(yL(d1,d2)
˜
yL(d1,d2) ·µd1·µd2,i,j )
where the filtered term is:
C
(yL(d1,d2)
˜
yL(d1,d2) ·µd1·µd2,i,j )
=∑
d1
∑
d2
pd1→i(−⃗ω)·pd2→j(−⃗ω)·yL(d1,d2)
˜
yL(d1,d2) ·µd1µd2. Now, let us consider the denominator of the back projection formula:
∑
L
ALV =∑
d1
∑
d2
AL(d1,d2),V ≈∑
i
∑
j
DL(i,j),V · C(µd1 ·µd2,i,j).
Summarizing, we can build the detector model into the back projection formula by computing a convolution for yy˜L(d1,d2)
L(d1,d2)·µd1·µd2 andµd1·µd2 before executing a geometric back projection.
The LOR filtering scheme is exactly the same as we perform in forward projection, just the direction vector needs to be reversed, this is why it is called Inverse LOR filtering.
6.3.2 Pre-computation
The input of our process is the crystal transport probability defined on the crystal structure, which has been computed by GATE in the following way [LCLC10]. Incident photons arriving from a direction of given inclination and azimuth angles at uniformly distributed points on the detector surface are simulated and the probabilities that this photon is absorbed in another crystal are computed. These probabilities can be visualized as a two dimensional image, where arrival probabilities are depicted by gray levels (see Figure 6.1 for an example). Also, this data can be interpreted as the (discretized) sub-critical probability density function of the offset
Figure 6.1: Flashing probabilities at right angle and in the case when both the inclination and azimuth angles are 40 degrees (logarithmic scale) [LCLC10].
vector between the photon impact point and the absorption point. There is a separate map for each discrete direction defined by a pair of inclination and azimuth angles.
The distinction according to azimuth angles is justified by the fact that crystals on the detector surface are squares and there are gaps filled by air between them.
The problem is that these images are too large to be sampled efficiently. So, during the pre-computation, we pre-generate relaxed MC sample sets that contain just a few samples, but their cumulative distribution is as close to the simulated distribution as possible. The relaxation step is included since we shall use these samples for high dimensional integration where the error is proportional to how well the density mimics the integrand and to the “distance” between the empirical cumulative distribution of the generated samples and the distribution used for the sample generation [SKS09].
As a result, we get the desired number of offset sample vectors for a given incoming inclination angle. The process is performed for all inclinations for which input measurement data exist, completing a sample offset set. For the actual reconstruction, several independent sample sets are generated for each inclination angle and sample number. By independent, we mean that the two sets are generated with independent random numbers.
6.3.3 Detector sampling during reconstruction
Offset sets were obtained by considering just one detector, so twoindependentsample sets should be used to sample the detector pair. According to the observations made in Section 2.2, we re-sample the system matrix in every iteration. We use two new independent offset sets for every iteration step, one for each of the two endpoints of LORs. Different LORs use the same offsets in an iteration step, which allows high performance on GPUs. If the iteration is longer than the half of the number of available sets, the sets are started to be used again.
6.3.4 Scatter compensation with detector response
So far we ignored scattering from the model when deriving the formulae for LOR filtering. If both scatter compensation and detector response are important (the measured object is big as in human PET and crystals are small as in small-animal PET), then both scattering and LOR filtering should be executed. However, we should be aware that it is only approximately feasible since these operators cannot be factorized. The explanation is that LOR space blurring L is also incident angle and energy dependent. For the direct component, we can use the direction of the LOR and the energy level of the electronϵ0 = 1. However, for the scattered component, the incident direction is the direction of the scattering points and not the other detector, and the energy level depends on where the annihilation happened and what the scattering angle was.
For example, for the case of single scatter simulation of Section 5.2 (S= 1), an approximate factoring would be the blurring of the LOR space once for each scattering point (their number
CHAPTER 6. DETECTOR MODEL 67 is typically a few hundred). When the line segments are combined (Step 4), then the number of hits on energy level 1 and the number of hits on lower level as well as the average energy are obtained (note that if annihilation happens on the line segment ending in this detector, then the energy level is 1; on the other hand, if annihilation happened in the other connected segment, then the energy level depends on the scattering angle as defined by the Klein-Nishina formula).
Having obtained these LOR images, a LOR filtering is executed for the two energy levels. Unlike filtering the direct contribution, the incident direction is computed from the direction between the detector and the scattering point. However, this is too costly computationally.
Thus, to allow the combination of expected hits due to direct contribution and scattering, we assume that the incident direction of scattering paths is similar to the LOR direction, and the incident energy is 1 relative to the energy of the electron if the energy is in the set energy window.