This chapter presented an efficient positron range compensation algorithm. The positron range calculation in heterogeneous material is decomposed to a series of positron range calculations in homogeneous materials, once for each material type of the examined object, and a final compositing step. Positron range in homogeneous material, in turn, is evaluated in the frequency domain applying Fast Fourier Transforms. The model runs on the GPU, providing positron range simulations at a negligible cost even for higher volume resolutions.
Chapter 4
Geometric projection
The majority of the measured coincident hits belong to direct photons that reach the detectors without scattering. For this direct component, the System Matrix (SM) is sparse and can be modeled by geometric projection, thus it is worth being treated independently of scattering. Effi-cient parallel implementation requires the geometric projection to be LOR driven in the forward projector and voxel driven in the back projector. Existing LOR driven forward projector meth-ods use varying sample number to evaluate line integrals and thus assign different computational load to parallel threads causing their divergence. Furthermore, they give analytic solutions to compute the direct contribution, which leads to a biased estimator and thus modifies the fixed point of the iteration. The LOR driven sampling scheme proposed in this chapter offers efficient parallel implementation using the same set of offsets in each thread. Furthermore, we re-sample the surfaces of the detectors in every iteration step of the ML-EM, and use a random offset for the line samples along the line to guarantee that every point that may correspond to a LOR is sampled with a positive probability.
The LOR driven approach may be wasting in the sense that it does not consider the emission density during sample generation. We propose a voxel driven geometric projection scheme that computes the contribution of a voxel to LORs, evenly sampling the detector surfaces. This allows the activity distribution to be taken into account in the forward projection, using importance sampling of the voxels. Furthermore, being a voxel centric approach, it provides an efficient parallel implementation of the back projector.
4.1 Previous work on geometric projection
4.1.1 Direct contribution between detectors
As we shall see later on in this chapter (see Equation 4.4), the direct contribution between two detectors can be expressed as a surface integral on the two detectors, where the integrand contains line integrals between the surface points:
˜ ygeomL(d
1,d2) =
∫
D1
∫
D2
⃗ z2
∫
⃗ z1
f(⃗l)dldz1dz2.
In general, these line integrals are computed as a sum of weighted samples taken by marching along the line, often referred to asray marching. Siddon’s algorithm [Sid85] is a usual choice for this task. Assuming a piece-wise constant finite-element approximation of the integrand f, the integral can be solved analytically as a sum∑
lVfV, wherelV is the length of the intersection of the ray with voxel V = (x, y, z) and fV is the voxel value inV. There are two major criticisms for this approach. The piece-wise constant approximation introduces an unrealistic discontinuity to the model [Jos82]. Furthermore, to evaluate the sum the intersections with the voxels are
37
computed along the ray which may need different number of loop cycles and divergent condi-tional instructions for different rays, thus, multiprocessor performance is degraded in a parallel implementation. Even so, there are existing GPU-based iterative PET reconstruction methods that apply Siddon’s approach [BVVC12]. Joseph’s method [Jos82] applies the trapezoidal rule to numerically estimate the integral and takes one sample per voxel along the ray, which are linearly interpolated from neighbouring values providing a smoother approximation of f. Ad-ditionally, the use of equidistant samples is more efficient on GPUs and the linear interpolation is provided with no additional cost, making Joseph’s approach more popular in GPU-based it-erative PET reconstruction [CM03, BS06]. However, different rays may intersect with different number of voxels causing varying number of loop cycles in threads. Thus, the sampling strategy we propose in Section 4.2.1 takes the same number of samples for every ray.
Surface integrals of the value obtained with the discussed line integral can be estimated by discrete line samples in a line driven method. For performance reasons, usually only a single line sample is used and thus the problem degrades to the computation of a line integral. The distance-driven approach [MB04] samples only one endpoint and simultaneously approximates the surface integral of the other endpoint and the line integral. Solid angle based methods[QLC+98] approxi-mate surface integrals. In general, the three integrals can be estiapproxi-mated analytically [MDB+08] or with Monte Carlo (MC) quadrature, or we can even mix the two approaches and some integrals are estimated with simple analytical formula while others are computed from random samples.
Integrating some variables analytically, we can increase the accuracy when low number of MC samples are used. However, analytical approximations have a deterministic error which makes the method biased, i.e. the error will not converge to zero when the number of MC samples goes to infinity. In order to get an unbiased estimator, Section 4.2.1 proposes random sampling of the detector surfaces and random shifting of the voxel samples along the lines, re-sampled in each iteration step of the ML-EM.
4.1.2 GPU-based projectors
An efficient, gather-style GPU implementation of the forward and back projectors must be LOR driven and voxel driven, respectively. In existing GPU implementations of a geometric forward projector [CM03, BS06, HEV+06], this principle is met. Despite the fact that it fits well to the massively parallel architecture, the LOR driven forward projection may be very inefficient since it completely neglects voxel values, i.e. in terms of importance sampling the sampling density does not mimic the emission density factor of the integrand. In the worst case, when the measured object is a point source, the algorithm wastes most of the samples traversing regions with zero activity while the sampling density around the point source is most likely to be insufficiently low. In Section 4.2.2, we propose a voxel driven projection that may consider the distribution of the activity. Being a voxel driven method, it would require an enormous amount of samples to accurately reconstruct large objects. Section 7.2 describes how to combine the benefits of the two sampling strategies, i.e. providing an accurate reconstruction for both point source-like and large objects with a reasonable amount of samples, according to multiple importance sampling.
Voxel driven back projection approaches [BS06, HEV+06, KY11a] search contributing LORs (for which the volume enclosed by the two endpoints of the LOR intersects with the voxel) by expressing LOR endpoints in polar coordinates, and sample them by looping through angles.
Since the axial cross sections of most scanners are not circles but polygons, this sampling of the detectors becomes uneven. The voxel driven projector presented in Section 4.2.2 samples the surface of the detector modules evenly and finds the other endpoint of the LOR via projection through the voxel, leading to more uniform sampling in LOR-space.
CHAPTER 4. GEOMETRIC PROJECTION 39