Evaluation of list-mode ordered subset expectation maximization image reconstruction for pixelated solid-state compton gamma camera with large number of channels

The Voxel Imaging PET (VIP) Pathfinder project intends to show the advantages of using pixelated solid-state technology for nuclear medicine applications. It proposes designs for Positron Emission Tomography (PET), Positron Emission Mammography (PEM) and Compton gamma camera detectors with a large number of signal channels (of the order of 106). For Compton camera, especially with a large number of readout channels, image reconstruction presents a big challenge. In this work, results are presented for the List-Mode Ordered Subset Expectation Maximization (LM-OSEM) image reconstruction algorithm on simulated data with the VIP Compton camera design. For the simulation, all realistic contributions to the spatial resolution are taken into account, including the Doppler broadening effect. The results show that even with a straightforward implementation of LM-OSEM, good images can be obtained for the proposed Compton camera design. Results are shown for various phantoms, including extended sources and with a distance between the field of view and the first detector plane equal to 100 mm which corresponds to a realistic nuclear medicine environment.


Introduction
This work intends to show the validity of the List-Mode Ordered Subset Expectation Maximization (LM-OSEM) algorithm with simulation data with a Compton gamma camera design as proposed by the Voxel Imaging PET (VIP) Pathfinder project [1]. Image reconstruction for a Compton camera is challenging because in order to estimate the radiation sources it has to combine Compton cones instead of lines-of-response (LOR) and the resolution of the cone estimation is not only depending on the energy and spatial resolution of the detector but, additionally and most importantly, also on the Doppler broadening effect [2] (explained below). The challenge is even bigger for the VIP Compton camera design because of its large number of independent readout channels. Compton gamma cameras use the kinematics of Compton scattering to localise the radioactive source. Gamma rays emitted by a radioactive source scatter inelastically in the scatter detector and, subsequently, are absorbed in the absorber detector (figure 1). The original gamma ray source is located on the surface of the Compton cone identified by the cone axis and apex, determined from the hit locations and the scattering angle θ , obtained from: where E γ is the known energy of the original gamma and E scatter is the energy deposited in the scatter detector. The Doppler broadening effect tells us that the electrons in the scattering material have an initial momentum, so that equation (1.1) is only approximately valid.

JINST 9 C04034
2 The VIP Compton camera The aim of the VIP project is to show that using pixelated solid-state technology with a large number of signal channels for nuclear medicine detectors like positron emission tomography (PET), positron emission mammography (PEM) and Compton cameras will improve both the energy and position resolution of the measurements with high efficiency compared to state-of-the-art crystal detectors [3,4]. For the readout a dedicated application specific integrated circuit (ASIC) is designed and produced by the project [5].
The basic building block of the VIP detector designs is the VIP unit module (figure 2) which is made of 4 CdTe pixelated detectors, of 200 voxels each where each voxel has a size of 1 × 1 × 2 mm 3 and is connected to its own independent readout channel for the energy, position, and arrival time. The unit module can be easily stacked to form large scale PET [6], PEM [7] or Compton camera applications. The scatter detector and the absorption detector of the VIP Compton camera are made of pixelated Si and CdTe sensors respectively (figure 1). The scatter detector has a parallelepiped shape of size 380 × 540 × 26 mm 3 and a nominal depth of interaction of 2 cm. The absorption detector of the Compton camera, at a distance of 10 cm from the scatter detector, has a parallelepiped shape of size 380 × 540 × 62 mm 3 and a nominal depth of interaction of 4 cm. The total number of independent readout channels amounts to more than 3 million. The expected efficiency of the Compton camera, for a 511 keV source at 100 mm distance from the scatterer, is about 6 cps/kBq with a signal purity (i.e., the ratio of true to total coincidences) of the order of 95% [8].
One of the basic characteristics of the VIP design is the flexibility in the number of basic modules that can be stacked to built the whole detector, depending on what is needed. After testing the ASIC this year, the VIP project will proceed with building and evaluating a PET prototype using a modified VIP unit module, and, subsequently, a simplified version of the VIP Compton camera design with a size of 10 x 10 cm 2 , as proofs of concept.

LM-OSEM
With Maximum Likelihood Expectation Maximization (MLEM), an image estimate is forward projected onto the detector and the projected data is compared with the real measured data. Next, a cost function is used to update the image estimate. These steps are repeated iteratively until a minimum error has been reached as depicted in figure 3. Ordered Subset Expectation Maximization (OSEM) [9] works the same as MLEM, but the algorithm is optimised by doing an update after each subset of the total amount of data has been processed. However, for a Compton gamma camera with a large number of channels, OSEM still requires too much memory and CPU-time consumption, because it uses a system matrix that maps probabilities for a signal coming from a particular voxel in the field of view (FOV) to be detected in a particular certain detector bin. Since the total number of detector bins equals all possible combinations of channels in the scatter detector and the absorber, the system matrix would be too large to be practical. With the data presented in list mode, LM-OSEM [10] sums over all detected events that have cones that intersect with each particular FOV bin instead of looping over detector bins. This significantly reduces the required CPU memory. Nevertheless, due to the large number of pixels in the FOV potentially lying on a particular cone, it still can be very time consuming.
The LM-OSEM iteration consists of two steps: • forward projection of an estimate of the image λ onto the detector: ∑ where t ik is the transition probability for event i to have originated from FOV bin k and M FOV i is the number of bins in the FOV that are intersected by the cone of event i.
• back-projection of the measured data, weighted with the forward projected data, providing an update correction for the image estimate.
These steps combine into the following equation: The transition probability t ik and the sensitivity s j (i.e., the probability for FOV bin j to produce a detected event) depend on physics considerations. The calculation of these variables would require the determination of the probability of each event to occur, starting from a particular FOV bin, by determining its particular cross-section. To reduce the cost in CPU time and memory, these variables are set to one for the results presented in this paper.
For the back-projection of the measured data, it is necessary for each event to find all FOV voxels that lie on its corresponding cone. For this, an accelerated back-projection algorithm is used [11], as illustrated in figure 4. Here, for each slice (with z = z s ), first the intersection of the cone with the edges of the FOV are calculated, using the following formula where (n x , n y , n z ) are the components of the unit vector along the Compton cone axis, (x 1 , y 1 , z 1 ) is the apex of the Compton cone and either x or y is known, depending on whether a intersection on a vertical edge or a horizontal edge is calculated: -3 - Subsequently, intersections with the two nearest adjacent pixel lines are searched for, and all pixels in between the previous solution and the next are marked as lying on the cone. The same procedure is repeated for each slice in the FOV.
Additionally, LM-OSEM has been implemented for the VIP PET and PEM detector design proposals. In these cases, since the algorithm deals with LORs instead of cones, a straight forward ray-box intersection algorithm is used [12] to find the FOV bins corresponding to each event. The evaluation of this implementation showed excellent results [13,14].

Results
All simulation data are obtained with the Geant4-based Architecture for Medicine-Oriented Simulations (GAMOS) package [15] and takes into account all contributions to the overall resolution, including the Doppler broadening effect. The simulation data is made, for all phantoms, by simulating a 511 keV gamma emitting radioactive source in an environment filled with air. The distance from the centre of the FOV to the front of the scatterer is 100 mm which is realistic for nuclear medicine purposes. The event selection is done as described in [8]. Point spread function (PSF). An 3D image defined by a FOV of 21 × 21 × 21 bins of size 0.5 mm 3 was reconstructed by our LM-OSEM implementation on simulated Compton camera data with a point-like source at position (1.5, 1.5, 0.0). The coincidence data set consisted of 7 million events and LM-OSEM ran for 11 iterations with the data set divided into 3 subsets. The 2D profiles of the final image along different orientations in figure 5 clearly show that the resolution along the z axis (i.e., perpendicular to the detector) is worse than along the x and y axes.
To obtain a better resolution along the z axis, we also simulated a Compton camera rotated 90 degrees around the y axis so that its planes are parallel to the zy plane, as illustrated in figure 6. In this case, the resolution in the yz is better than along the perpendicular x axis.
Another thing to notice in figure 5 are the artificial signal in the edges and, especially, the corners of the images. These are due to background events which, although they are present in the entire 3D FOV, for cones which only intersect with the FOV corners will get trapped in the only FOV pixels available for them. For the results presented in all following figures, the outer edges of the images are ignored.
-4 -  Figure 7 shows the line profiles measured with the Compton camera in two different orientations. The FWHM of the PSF is 1.9 mm in the x direction and 2.9 mm in the y direction, with the Compton camera oriented along the xy plane. With the Compton camera rotated 90 degrees around the y axis and oriented along the yz plane, the PSF (FWHM) is 1.9 mm in the z direction and 2.9 mm in the y direction. The bigger spatial resolution in the y direction reflects the detector voxels size, which is 2 mm in the y direction and 1 mm in the x and z direction. 3D Cube phantom Figure 8 shows a 3D cube phantom, consisting of 8 spheres on the corners of a cube. The spheres have a diameter of 2 mm and a centre-to-centre distance of 5 mm. Figure 9 shows the 2D images and profiles along the x and y axes for slices at z = 2.5 mm and -2.5 mm. As explained before, to get a good resolution in the z direction, the Compton camera should be rotated 90 degrees around y. Figure 10 shows the 2D images and profiles along the y and z axes for slices at x = 2.5 mm and -2.5 mm, taken with the rotated Compton camera.
Horseshoe phantom. An important requirement for nuclear medicine applications is the ability to correctly reproduce extended source objects. For this purpose, a horseshoe shaped phantom is used, as depicted on the left in figure 11. The width and the thickness of the horseshoe are both 1 mm and the gap size is 4.5 mm.
The image on the right in figure 11 shows very good results can be obtained with LM-OSEM on Compton camera data with the horseshoe phantom. These results were obtained with a data set of 20 million events and LM-OSEM ran with 5 subsets, for 8 iterations, with a FOV of 40 × 40 × 1 bins of size 0.5 × 0.5 × 3 mm 3 . The result in figure 11 shows a "hot spot" on the left hand side of the phantom. This is caused by the detector voxels size, which is bigger for y (2 mm) than for x (1 mm), causing a bigger spread in the image in y than in x. When the voxels would have a size of 1 × 1 × 1 mm 3 this effect would disappear.
-5 -   Derenzo phantom Another phantom with extended source objects is the Derenzo phantom, which serves as a measure of the spatial resolution of the detector. The Derenzo phantom consists of 5 segments, each containing rods with equal activity, of length 12 mm and with varying diameters and distances as indicated in figure 12.
-6 -  Figure 13 shows the good results obtained with LM-OSEM on Compton camera data with this phantom. The rods with diameter 1.5 mm and 6 mm distance between the rod centres are distinctly reconstructed, as also can be seen from the line profile on the right of figure 13. These results were obtained with a data set of 70 million events, and LM-OSEM ran with 70 subsets, for 2 iterations, with a FOV of 160 × 160 × 1 bins of size 0.5 mm × 0.5 mm × 12 mm.

Conclusion
The performance of a straigtforward implementation of LM-OSEM (i.e., without taking into account the transition and absorption probabilities of the Compton scattering process) has been studied on simulation data with the VIP Compton camera design and various phantoms. The simulations include the Doppler broadening effect and the distance between the FOV centre and scatterer was kept at a distance of 100 mm, which is realistic for nuclear medicine purposes.
For the PSF with a point-like source, LM-OSEM achieves values of the order of 2 mm (FWHM). By combining the data taken with the Compton camera positioned in two orthogonal orientations a good spatial resolution can be obtained in all directions. Similarly, LM-OSEM can produce 3D images of spatially distributed 511 keV gamma point-like sources with a distance of 5 mm. Also with more complicated extended sources excellent results were obtained. Results with a horseshoe shaped phantom showed that both the shape and the 4.5 mm gap can be very well reconstructed. With a Derenzo phantom rods with a diameter of 1.5 mm and a centre-to-centre distance of 6 mm were succesfully reconstructed.
Future work will include the improvement of the spatial resolution in all directions by further studies on rotating the Compton camera, the study on how to approximate the transition probability t ik and the sensitivity s j of the LM-OSEM algorithm and a systematic comparison of the performance of LM-OSEM and other algorithms for Compton camera (e.g. OE, as reported in [17]) with various phantoms. Additionally, it has to be studied how to diminish the artefacts in the edges and, especially, the corners of the images.