3D Compton image reconstruction method for whole gamma imaging

Compton imaging or Compton camera imaging has been studied well, but its advantages in nuclear medicine and molecular imaging have not been demonstrated yet. Therefore, the aim of this work was to compare Compton imaging with positron emission tomography (PET) by using the same imaging platform of whole gamma imaging (WGI). WGI is a concept that combines PET with Compton imaging by inserting a scatterer ring into a PET ring. This concept utilizes diverse types of gamma rays for 3D tomographic imaging. In this paper, we remodeled our previous WGI prototype for small animal imaging, and we developed an image reconstruction method based on a list-mode ordered subset expectation maximization algorithm incorporating detector response function modeling, random correction and normalization (sensitivity correction) for either PET and Compton imaging. To the best of our knowledge, this is the world’s first realization of a full-ring Compton imaging system. We selected 89Zr as an imaging target because a 89Zr nuclide emits a 909 keV single-gamma ray as well as a positron, and we can directly compare Compton imaging of 909 keV photons with PET, a well-established modality. We measured a cylindrical phantom and a small rod phantom filled with 89Zr solutions of 10.3 MBq and 10.2 MBq activity, respectively, for 1 h each. The uniform radioactivity distribution of the cylindrical phantom was reconstructed with normalization in both PET and Compton imaging. Coefficients of variation for region-of-interest values were 4.2% for Compton imaging and 3.3% for PET; the difference might be explained by the difference in the detected count number. The small rod phantom experiment showed that the WGI Compton imaging had spatial resolution better than 3.0 mm at the peripheral region although the center region had lower resolution. PET resolved 2.2 mm rods clearly at any location. We measured a mouse for 1 h, 1 d after injection of 9.8 MBq 89Zr oxalate. The 89Zr assimilated in the mouse bony structures was clearly depicted, and Compton imaging results agreed well with PET images, especially for the region inside the scatterer ring. In conclusion, we demonstrated the performance of WGI using the developed Compton image reconstruction method. We realized Compton imaging with a quality approaching that of PET, which is supporting a future expectation that Compton imaging outperforms PET.


Introduction
Compton imaging or Compton camera imaging is a method to image gamma ray sources by measuring a Compton scatter and a photoelectric absorption for each decay in terms of positions and deposited energies. The source position for a measured event is identified on a surface of a cone determined by the Compton kinematics with the measured information. For a point source, the source position is estimated by drawing and overlaying multiple cones. For more complicated source distributions, analytical and iterative image reconstruction methods were developed (Singh and Doria 1981, Singh 1983, Hebert et al 1990, Cree and Bones 1994, 2000, Parra 2000, Tomitani and Hirasawa 2002.
The Compton imaging method was originally proposed for solar neutron imaging (Pinkau 1966, White 1968). The first practical system was developed as a Compton telescope to measure radiation from outer space (Schönfelder et al 1973(Schönfelder et al , 1984(Schönfelder et al , 1993. Since then, several nondestructive assay and environmental radiation measurement systems have been developed (King et al 1994, Mihailescu et al 2007, Jiang et al 2016, Muraishi et al 2020. The first clinical application was proposed as a gamma camera for nuclear medicine imaging (Todd et al 1974, Everett et al 1977. Toward practical use of the Compton imaging for nuclear medicine, various studies of theoretical analysis, instrumentation development, and reconstruction algorithm development were conducted (Singh 1983, Singh and Doria 1983, Singh and Brechner 1990, Martin et al 1993, Gormley et al 1997, Leblanc et al 1998, 1999, 2000, and such work was continued. Several researchers applied Compton imaging for biomedical applications with living subjects (Motomura et al 2007, 2013, Kabuki et al 2010, Takeda et al 2012, Kishimoto et al 2017, Sakai et al 2018, Nakano et al 2019. Most previous studies were single view imaging with a single detector or 3D imaging by rotating one or two detectors, which might have limited the projection views and the radiation counts in image reconstruction. To the best of our knowledge, nobody has developed a full-ring Compton imaging system. In nuclear medicine and molecular imaging, Compton imaging of various energy gamma rays has been demonstrated but imaging quality has been far from that of typical positron emission tomography (PET). In contrast to PET, which can only visualize positron emitters, Compton imaging has a potential to visualize a variety of gamma ray sources. However, the advantages of Compton imaging in nuclear medicine and molecular imaging have not been demonstrated yet. Therefore, the aim of this study was to compare Compton imaging with PET by using the same imaging platform of whole gamma imaging (WGI). WGI is a concept that combines PET with Compton imaging by inserting a scatterer ring into a PET ring (Yamaya et al 2017, Yoshida et al 2017a. This concept utilizes diverse types of gamma rays for 3D tomographic imaging by measuring single-gamma rays with the Compton imaging technique in addition to PET coincidence measurement. A WGI system works as both a PET system (PET mode) and a Compton imaging system (Compton imaging mode). Similar double ring designs have been proposed and studied for Compton imaging applications (Bolozdynya et al 1997, Liang and Jaszczak 1990, Valentine et al 1997, Rogers et al 2004 and for high-resolution PET systems (Park et al 2007a, 2007b, Yin et al 2014, but hybrid PET and Compton imaging systems have not been studied yet except by Yoshino 2017 andShimazoe et al 2020. In this study, we realize the world's first full-ring Compton imaging based on the WGI geometry. For 3D tomographic imaging, image reconstruction with appropriate corrections is important for achieving high quality images. Therefore, we develop an image reconstruction method incorporating detector response function (DRF) modeling, random correction and normalization for the WGI Compton imaging. After confirming the developed Compton imaging method with phantoms, the WGI prototype is applied to the first animal imaging. For this purpose, the diameter of the scatterer ring is half-sized to better fit mouse imaging as spatial resolution of Compton imaging is basically degraded in proportion to the source-detector distance. For a positron emitter, we use 89 Zr because 89 Zr emits a 909 keV gamma ray after a positron decay. We think that 89 Zr should be an optimum isotope to compare Compton imaging with PET. The half-life of 89 Zr is 3.27 d, and it decays to 89m Y via positron emission (23%) and via electron capture (77%). Then, 89m Y decays to a stable isotope of 89 Y via single-gamma emission with a half-life of 15.7 s. Therefore, we can image the same distribution independently by either Compton imaging or PET.

Modification of a WGI prototype for small animal Compton imaging 2.1.1. Gantry design and detector arrangement
We modified the WGI prototype (Yoshida et al 2017a(Yoshida et al , 2020 for the small animal experiment by half-sizing the diameter of the scatterer detector ring. The prototype consisted of a PET ring (absorber ring) with the diameter of 660 mm and the scatterer ring insert with the diameter of 94 mm (figure 1). The PET ring was a whole-body depth-of-interaction (DOI) PET scanner (Akamatsu et al 2019). The axial lengths of absorber and scatterer rings were 214 mm and 52 mm, respectively. The absorber ring had four rings of 40 4-layer DOI detectors each consisting of an array of 16 × 16 × 4 Zr-doped gadolinium oxyorthosilicate (GSOZ) (Shimura et al 2006) crystals (2.8 × 2.8 × 7.5 mm 3 ) and a 64-ch photomultiplier tube (R10552-100-M64) (Hirano et al 2014). The scatterer ring had two rings of ten scatterer detectors, each consisting of an array of 24 × 24 gadolinium aluminum gallium garnet (GAGG) (Kamada et al 2012) crystals (0.9 × 0.9 × 6.0 mm 3 ) Figure 1. Design sketches of the modified WGI prototype shown as a side view (a) and a front view (b); photographs of the scatterer detector and the absorber detector (c), and the developed prototype (d).
and an array of 8 × 8 multi pixel photon counters (MPPCs) (S14161-3050HS-8, Hamamatsu Photonics K. K.) with the pixel size of 3.0 × 3.0 mm 3 . The design of the scatterer detectors strove toward high spatial resolution and high energy resolution with the high light output of GAGG. To discriminate 0.9 mm crystals by 3.0 mm MPPC pixels, a 2.5 mm thick light guide was inserted between the crystal array and the MPPC array. In order to compensate for saturation of MPPC signals, energy information for the scatterer detectors was obtained by applying a saturation correction method which used measured energy spectra of 133 Ba (288 keV), 68 Ge (511 keV) and 137 Cs (662 keV) sources. We should note that, unfortunately, the energy resolution was found to be worse than we expected due to an unoptimized pulse integration time in the data acquisition (DAQ) system.

Data acquisition and event selection for 89 Zr experiment
The DAQ system of the WGI prototype recorded all single events during measurements. Each single event contained a detector index, crystal position information, 8-bit energy information, and a time stamp of 500 ps resolution. We applied an event selection process such as energy discrimination and coincidence detection as a post postprocessing with offline software. Figure 2 shows a flow chart for the event selection process for PET and Compton imaging measured for 89 Zr. As the energy resolutions of the scatterer detector and the absorber detector were 17% and 14%, respectively, at 511 keV (Yoshida et al 2017a(Yoshida et al , 2017b(Yoshida et al , 2020, we set the energy windows for PET events to 350-650 keV for the scatterer ring and 400-600 keV for the absorber ring. After energy discrimination, software coincidence was done with a relatively wide timing window of 50 ns, because no precise timing correction was applied. We should note that there are three types of PET events: (1) scatterer and scatterer, (2) scatterer and absorber, and (3) absorber and absorber. All these events were used for PET image reconstruction. Inserting a detector ring of smaller crystals may improve spatial resolution as well (Park et al 2007a, 2007b. For Compton imaging of 909 keV photons, the energy windows were set to 50-350 keV for the scatterer ring and 550-850 keV for the absorber ring. The coincidence events between the scatterer ring and the absorber ring were extracted with the timing window of 50 ns, and another energy window of 800-1000 keV was applied for the total energy to extract 909 keV photons. The PET and Compton events were separately extracted from the same set of singles events. Delayed coincidence events for both modalities were also extracted by shifting the coincidence window 256 ns.

3D Compton image reconstruction method 2.2.1. Detector response model of Compton imaging
We implemented the list-mode ordered subset expectation maximization (OSEM) algorithm (Shepp and Vardi 1982, Hudson and Larkin 1994, Reader et al 1998, Parra and Barrett 1998, Rahmim et al 2005 for either Compton imaging and PET. One of the key components in implementation is a DRF incorporating the angular blurring effect. We should note that, for high-quality image reconstruction, it is necessary to model the DRF representing the actual system response as accurately as possible. We modeled the DRF for Compton imaging as follows.
The scatter angle θ is calculated using the Klein-Nishina equation as

PET mode Compton imaging mode (909 keV)
Compton events PET events Apply energy window for total energy 800 keV < Es + Ea < 1000 keV where E γ is the original gamma ray energy emitted from radioisotopes, E s is the energy deposited at the scatterer, m e is the electron mass, and c is the speed of light. The measured scatterer energy contains uncertainty due to the energy resolution. The energy resolution depends on the deposited energy. For the DRF modeling, we assumed that the energy deviation follows a Gaussian distribution, with the full width at half maximum (FWHM) modeled as follows: where α and β are parameters determined by fitting this model to measured FWHM values for several deposited energies. Figure 3 shows the DRF model. We modeled the blurring on the Compton cone due to the angular deviation as an asymmetric Gaussian function with standard deviations directly converted from increased and decreased energy values with the half width at half maximum (HWHM) %HWHM(E s ) = %FWHM(E s )/2. The angle with the increased energy value is calculated as ) m e c 2 + 1 ) , and that with the decreased energy value is calculated as ) m e c 2 + 1 ) .
We let x be a position in the imaging domain (e.g. a voxel center), s be the detection position at the scatterer, and a be the detection position at the absorber. Then, the standard deviations corresponding to the energy increase and decrease notated as σ + and σ -, respectively, are calculated as follows: where n is the directional vector from scatterer position s to voxel position x, θ x is the angle between n and the line from s toward x, x r is the distance from the cone surface to x, x c is the distance between s and x projected on the cone surface. The blurring kernel is represented as Finally, the DRF of list-mode data consisting of s, a, E s for the position x is represented with incorporation of distance weights as follows: The DRF model does not include position uncertainty of s and a determined as fixed points in the scatterer and absorber detectors. We assumed that the effect of the position uncertainty is small compared with the angular blurring due to the energy deviation, but the DRF can be extended to include the position uncertainty effect by placing multiple sampling points in crystals and drawing multiple cones for a list-mode event at the cost of an increased computational burden. Figure 4 shows Compton cones drawn by the DRF model for choosing several cone angels for 909 keV and 511 keV. Based on the energy resolution of the WGI prototype, α and β were set to 900 and 15, respectively.   For each image pair, the upper shows the cross section perpendicular to the direction from the absorber detection point to the scatterer detection point, and the lower shows the cross section parallel to the direction. The scatterer detection point is at 50 mm below (the radius of the scatterer ring, 47 mm + half size of the crystal length, 6 mm/2) the center of bottom image in each pair and the absorber detection point is on the extended line from the center to the scatterer. Each upper image also shows the energy deposited at the scatterer (Es), the energy deviation (%FWHM(Es)), and perturbed angles by energy decrease (θ − θ − ) and increase (θ + − θ).

List-mode Compton image reconstruction algorithm
To formulate a reconstruction algorithm, we should define a system matrix relating measurement domain and imaging domain. We let i be an index of a measured data bin consisting of a set of discretized parameters s, a, E s and j be a voxel index positioned at x. The system matrix element a ij incorporating DRF can be calculated with equation (12).
Correction for the random coincidence events (random correction) is done by the delayed-coincidence method, which is the most common technique for PET scanners (Rahmim et al 2005). The delayed-coincidence events (delays) are acquired by setting a delayed-coincidence window having the same width as the coincidence window for extracting prompt events (prompts) containing true, object scatter, and random events.
The implemented list-mode OSEM algorithm with the random correction and the normalization, which is a method to correct the sensitivity difference between the calculation model and the actual system, is summarized as follows: where f is the reconstructed image; C is the global sensitivity image; the superscripts k and l indicate the main and sub-iterations, respectively; the subscripts i, j, and t indicate measured data bin, voxel, and coincidence event, respectively; J and L are the numbers of voxels and subsets, respectively; S l is the set of list-mode data in the lth subset; and We note that equation (14) is equivalent to a histogram-based algorithm subtracting delays from projection data bins for each subset (Rahmim et al 2004). For normalization, we assumed sensitivity uncertainty is due to characteristics of the detector elements and only diagonal sensitivity matrix factors are required rather than all matrix elements for element-wise product. Therefore, normalization can be summarized as the sensitivity image C. Normally, back-projection of all data bins with normalization factors is required. Although the components of the list-mode data are discretized through a DAQ system for computer implementation, the total number of data bins is huge as it is a multiplication of the bin numbers of scatter position (number of scatterer crystals), absorption position (number of absorber crystals), and energy. Especially, it is difficult to arrange the energy information into a histogram bin because each detector element has a different correction factor and resampling causes energy resolution degradation. Therefore, we used the list-mode data directly even for the normalization to calculate the global sensitivity image: where p is a known radioactivity distribution of a normalization phantom, and M is the set of list-mode data acquired by measuring the normalization phantom. For PET scanners, it is common to use line sources rotating around the field of view (FOV) having a cylindrical shape. The time integrated activity distribution of the rotating line sources is equivalent to the shape of a hollow cylinder (Mizuta et al 2010). For the full-ring Compton imaging, we proposed to use a normalization phantom with a hollow cylinder shape (hollow cylindrical phantom) filled with the same radioisotope as injected for subjects and with the radioactivity distribution surrounding the FOV. The measured data are compared with the forward projection with the system matrix and back-projected onto the image domain. We should note that this is equivalent to the direct normalization method used for PET normalization. At this stage, we have not yet implemented attenuation and scatter corrections because these effects are negligible for small subjects with high energy gamma rays such as those emitted from 89 Zr.

Experiment
We conducted phantom and small animal experiments of the WGI prototype with 89 Zr. Because a 89 Zr nuclide emits both a positron and a single-gamma ray of 909 keV, the same radioactivity distribution is obtained for both PET and Compton imaging and we can directly compare PET and Compton imaging of WGI.

Cylindrical phantom experiment
We measured a cylindrical phantom for the uniformity evaluation. The container of the phantom was made of low-density polyethylene with a wall thickness of 0.5 mm. The inner diameter was 38 mm and the axial length was 60 mm. We filled the phantom with a 89 Zr solution with the total radioactivity of 10.2 MBq. The phantom was placed at the center of the FOV aligned to match the cylinder axis with the axial direction of the prototype. We measured the phantom for 1 h. To confirm the effectiveness of the normalization method, we reconstructed images with and without normalization. In the case without normalization, a fixed value for each measurement bin was back-projected onto the image domain. As the number of measurement bins for the Compton imaging was huge, we coarsely sampled for the energy bin and the absorber detector position. Image uniformity was evaluated for reconstructed images by calculating a coefficient of variation (COV) for 63 cylinder-shaped regions of interest (ROIs) each with the diameter of 10 mm and the axial length of 5 mm placed at nine axial locations (centering at z = −20 mm to 20 mm with 5 mm steps) and seven transaxial locations (one at the center and the others surrounding it with a 1 mm separation interval). The COV was calculated as a percentage of the standard variation divided by the mean of the summed voxel value for each ROI.

Small rod phantom experiment
We measured a small rod phantom for spatial resolution demonstration in iterative image reconstruction.
The phantom consisted of an outer hollow cylinder and an inner solid cylinder with a diameter of 36.1 mm and a length of 17.8 mm. The inner solid cylinder had cylindrical hole clusters that could be filled with a radioactive solution. The hole diameters were 1 mm, 1.6 mm, 2.2 mm, 3.0 mm, 4.0 mm, and 4.8 mm, and the holes in each cluster were separated by a distance equal to the hole diameter. The phantom was placed at the center of the FOV with the rod direction aligned along the axial direction. The phantom was filled with a 89 Zr solution with total radioactivity of 10.3 MBq. We measured the phantom for 1 h. Separation of rods in reconstructed images was evaluated.

Small animal experiment
For the small animal experiment, we injected 9.8 MBq 89 Zr oxalate into an ICR mouse (male, 5 weeks old, Japan SLC) 23 h before measurement. We conducted a 1 h long measurement, in which the mouse torso was positioned inside the scatterer ring, while the head and tail parts were located outside the scatterer ring as the mouse was longer than the axial length of the scatterer ring. The FOV larger than the detector coverage can be obtained in Compton imaging. We should note that the axial length of the absorber ring (PET ring) was long enough to cover the entire mouse body. The bed height was set at an elevated position to place the mouse closer to the scatterer ring because the spatial resolution of the Compton imaging is better when the subject is closer to the scatterer detector. In addition, we measured another 1 h only for Compton imaging by placing the upper half of the mouse body inside the scatterer ring to compare characteristics of Compton imaging inside and outside the scatterer ring. The radioactivity distributions in Compton images were compared with a PET image, and noise in images were visually evaluated. This experiment was approved by the local ethical committee of the National institute of Radiological Sciences, the National Institutes for Quantum and Radiological Science and Technology (Approved number: 07-1064-24) and carried out in compliance with the institute's animal experimentation guidelines.

Normalization measurement
For normalization, we measured a hollow cylindrical phantom. The phantom was designed to have radioactivity only on the surface of a cylinder with no radioactivity nor attenuation material inside. The imaging region with normalization would be inside the hollow cylindrical phantom. The inner diameter and axial length of the annulus were 80 mm and 274 mm, respectively, and the thickness of the space to be filled with radioactive solution was 3 mm. The phantom was made of polymethyl methacrylate with a thickness of 2 mm for the curved walls. The phantom was filled with 89 Zr solution with the radioactivity of 5.3 MBq. The phantom was inserted into the scatterer ring to be placed at the center of the gantry. We measured the phantom for 24 h to acquire enough counts for the normalization method. The 89 Zr distribution with the hollow cylinder shape was assigned to the known radioactivity p in equation (15).

PET and Compton image reconstruction
For each experiment, PET events and Compton events were extracted by offline processing with the coincidence time window of 50 ns and the delayed coincidence window shifted 256 ns, and images were reconstructed by the list-mode OSEM. The numbers of subsets and iterations were respectively 16 and 10 for PET mode and 16 and 50 for Compton imaging mode. For the PET image reconstruction, DRF for each list-mode event was modeled by a simple Gaussian function defined on the distance between the voxel center position and the LOR with the FWHM determined by considering the average size of the crystals composing the LOR. For the Compton image reconstruction, parameters α and β in equation (2) modeling energy resolution were set to 900 and 15, respectively, by taking the energy resolution of the scatterer detectors measured as 17% at 511 keV (Yoshida et al 2020) into account. Attenuation and attenuation corrections were not applied for the two reconstruction methods. Table 1 shows true + scatter (object scatter) counts estimated by subtracting the delayed coincidence event counts from the prompt coincidence event counts for each measurement. For PET, coincidence events within the scatterer and scatterer (S-S), those between the scatterer and the absorber (S-A), and those within the absorber ring (A-A) were extracted. Because the scatterer ring had thinner crystals and a shorter total length in axial direction than the absorber ring, the number of coincidence events within the scatter ring (S-S) was only 2% of that within the absorber ring (A-A). We should note that we used all of the event types together in the PET image reconstruction. The number of Compton events was 10%-15% of the PET events, although the number of single-gamma ray emissions for 89 Zr was 4 times larger than the number of positron emissions. We should also note that events of 511 keV photon pairs in that one side or both sides were scattered at the scatterer ring and absorbed at the absorber were not extracted in the current scheme. The random fractions were 14%-20% for the PET events and 30%-40% for the Compton events. Figure 5 shows the reconstructed images for the cylindrical phantom experiment. To show the effect of the normalization method, we did reconstruction with and without normalization. Without normalization, the PET image showed ring artifacts in the transaxial slices and stripe patterns in the sagittal slice. On the other hand, the Compton image did not show any specific pattern of artifacts but the axially center region had  higher values than the off-center regions. This was because the DRF for Compton imaging was much more blurred and widely distributed than that for PET. The blurred and widely spread DRF smoothed out the small structures in the sensitivity distribution. With normalization, both PET and Compton images were successfully reconstructed with uniform radioactivity intensity throughout the phantom. The uniformities of the radioactivity intensity for ROIs calculated as the COV without and without normalization were respectively 10.6% and 3.7% for the PET images and respectively 8.6% and 4.8% for the Compton image. For the Compton imaging, the radioactivity distribution outside the scatterer ring was shrunk in the axial direction compared with the PET image, because only obliquely limited information was available for the region outside the scatterer ring. The nozzle part of the cylindrical phantom shown in the PET image was shortened in the Compton imaging. Figure 6 shows imaging results of the small rod phantom experiment. The PET image resolved 2.2 mm rods clearly. The Compton image resolved 3.0 mm rods for the peripheral region, which was also confirmed in the profile. For the center region, the rods were not resolved clearly although the radioactivity distribution agreed with the PET mode results. We concluded that the region closer to the scatterer detector had better spatial resolution for Compton imaging. Figure 7 shows the reconstructed images for the small animal experiment. Because the absorber ring was axially long enough relative to the mouse length, the entire mouse body was clearly imaged with uniform quality by PET. The PET image showed that 89 Zr was assimilated in the mouse bones, and the distribution was clearly depicted. In the Compton mode, on the other hand, even though the scatterer ring with the axial length of 52 mm was shorter than the mouse length, the distribution outside the ring was also reconstructed and almost the entire body was imaged. On visual inspection, the overall contours were almost the same between the PET and Compton images. However, as imaging conditions were different, we observed different image qualities between regions inside and outside the scatterer ring. The image for the inside region agreed well with the PET image. On the other hand, we observed high accumulation on the top of the head located in the outside region, which was not observed in the PET image. Figure 8 shows the Compton imaging result when the upper half of the mouse body was placed inside the scatterer ring. The high accumulation observed when located outside ( figure 7(b)) was not observed when located inside (figure 8). We noticed that the two high accumulations on each leg were clearly depicted in figure 7(b) while they were combined in figure 8 with blurred image quality. Even though the image quality inside the scatterer ring in the Compton image was better than that outside, the image quality was relatively lower than for the PET image; for example, focusing on the mouse skull bone, we saw that it was difficult to identify the orbit (cavity of the eyes) and the mandible (inferior jaw bone) on the Compton image while we could clearly identify them on the PET image. We should note that, although the image quality was not as good as the PET image, the WGI Compton imaging could provide highly accurate measurements for the region inside the scatterer ring.

Discussion
In this work, we compared Compton imaging with PET by using the same imaging platform of WGI. We selected 89 Zr as an imaging target because a 89 Zr nuclide emits a 909 keV single-gamma ray as well as a positron, and we could directly compare Compton imaging of 909 keV photons with PET, which is a well-established modality. As shown in figure 6, the small rod phantom experiment showed that the WGI Compton imaging had spatial resolution better than 3.0 mm at the peripheral region although the center region had lower resolution. PET resolved 2.2 mm rods clearly at any location. However, in the mouse measurement with 89 Zr oxalate shown in figure 7, Compton imaging results agreed well with PET images. The phantom and the small animal experiments with 89 Zr demonstrated the high quality Compton imaging performance of the WGI prototype with the developed image reconstruction method incorporating DRF modeling, random correction and normalization methods. We think there were four key points in our successfully getting Compton images with quality approaching that of PET images. The first key point was the high energy of the gamma ray emitted from 89 Zr. It is well known that the higher the energy the better the angular resolution. Spatial resolution of Compton imaging rapidly degrades when the gamma ray energy is lower than around 500 keV (Bandstra et al 2006, Kurosawa et al 2009, Takeda 2009). Our resolution model also indicated that the DRF for 909 keV gamma rays was significantly sharper than that for the 511 keV ( figure 4). In addition, for high energy gamma rays, the effect of Doppler broadening can be significantly reduced compared with the cases of low energy gamma rays (Zoglauer and Kanbach 2003). Other candidate radionuclides emitting high energy gamma rays suitable for Compton imaging would be 95 Tc (gamma ray of 766 keV and half-life of 20 h) and 96 Tc (gamma rays of 778, 812, and 850 keV and half-life of 4.3 d), which can substitute a standard SPECT nuclide, 99m Tc (Hayakawa et al 2018).
Secondly, the DRF model used for image reconstruction was also important to achieve high spatial resolution. In this study, we derived the DRF model of Compton imaging by assuming the blurring effect depends only on energy uncertainty. However, the results inferred that the derived model approximated the actual system with sufficient accuracy.
Thirdly, normalization was essential to maintain the uniformity in images. It is often difficult to model an imaging system accurately, and the detection probability on each system matrix element for Compton imaging can significantly vary due to many factors such as characteristics of crystals, the quantum efficiency of photosensors, the incident angle and the scatter direction. Normalization is necessary to correct such variation, but appropriate normalization methods have not been investigated for Compton imaging. In this study, therefore, we transferred a method used in PET image reconstruction into Compton imaging. Although there may be several choices for an appropriate normalization phantom for Compton imaging, we used the same phantom for both modes for simplicity. In PET, typical iterative algorithms make a histogram of measured data and calculate normalization factors for all the histogram bins to be back-projected for the calculation of the global sensitivity image. However, it is difficult to make a histogram for a Compton imaging system because the size of the histogram, which is the multiplication of the numbers of scatterer detector elements, absorber detector elements and energy bins, becomes huge. Therefore, we developed a method to calculate the global sensitivity image directly from the list-mode data measured for normalization. We note that, for PET image reconstruction, random sampling methods of the histogram bins to calculate the global sensitivity image were proposed and proven to be effective for normalization (Carson et al 2003, Qi 2006. The list-mode acquisition is equivalent to the random sampling of the histogram bins according to the detection probabilities. In this study, we experimentally demonstrated the effectiveness of the normalization method directly using the list-mode data. As shown in figure 5, the uniformity evaluated by the COV for ROI values was 3.7% for PET and 4.8% for Compton imaging, which indicated that the implemented correction methods were effective for achieving highly accurate Compton images comparable to PET images. The last key point was the full-ring geometry, which enabled enough sampling of projection view angles as well as increasing the number of counts. In the small animal study, almost the entire body of the mouse was imaged by Compton imaging even though the size of the mouse was longer than the scatter ring width and some part of the body was located outside the ring. However, we observed degraded image quality for the region outside the scatterer ring. Also, in the case of the cylindrical phantom study, the region seemed to be shrunk as the nozzle part was shortened. The difference in imaging conditions for the region inside and outside the scatterer ring can be discussed with the data sufficiency condition as follows. Assuming that the line integrals starting from the scatterer detector elements can be reconstructed (Smith 2011) with sufficient accuracy with the WGI geometry, the position of the scatterer detector elements corresponds to a source trajectory for the cone-beam reconstruction using a circular orbit. Then, projection data are substantially complete for the region surrounded by the trajectory (inside the scatterer ring) because the region satisfies the data sufficiency condition that any plane passing through this part of the object intersects the source trajectory at least once (Tuy 1983, Smith 1985, 1990, Tang et al 2018. The projection data are incomplete for the region not surrounded (outside the scatterer ring) because the region does not satisfy the condition; therefore, the reconstructed image showed image quality degradation and shrinkage for this region. These considerations indicated that covering the target region of the subject by the scatterer ring is important for high-precision imaging, and it is desirable to develop an axially long scatterer ring for covering the subject.
Our next work will focus on improving Compton imaging performance so as to outperform PET. The energy resolution of the scatterer detector of the WGI prototype was 17% at 511 keV, which was not as good as typical Compton imaging systems. Nevertheless, we could achieve spatial resolution that was better than 3 mm. If we could achieve state-of-the-art energy resolution for the scatterer ring, the spatial resolution would be much better. In addition, the sensitivity of Compton imaging also needed improvement. The number of Compton events was only 10%-15% of the PET events (table 1), even though the number of single-gamma ray emissions was four times larger than that of the positron emissions (branching ratio of 23%). This was because the scatterer ring was shorter than the imaging subjects and the thickness of the crystals for the scatterer as well as the absorber was not optimized for 90 keV gamma rays in terms of the sensitivity. Therefore, we will develop a scatterer ring having high energy resolution with an extended axial length to cover the entire small animal body with the goals of high spatial resolution and high sensitivity as well as uniformity improvement. Our future development and investigation of improved scatterer detectors will include GAGG-based detectors with a DAQ system optimized the pulse integration time for high energy resolution and those with DOI measurement capability for high spatial resolution. Also, the use of other scatterer materials such as Si, CdTe, and CZT (Krimmer et al 2015, 2018, Turecek et al 2018, Draeger et al 2018 will be investigated.
Improvement of the energy resolution may enable the use of typical SPECT tracers for Compton imaging. In this study, we only applied the WGI Compton imaging for 909 keV gamma rays of 89 Zr because the current energy calibration setting targeting high energy gamma rays was not suitable for lower energy gamma rays. Although the Compton imaging is suitable for higher energy radionuclides, it can be expected that lower energy radionuclides such as SPECT tracers will play an important role in nuclear medicine, and there will be a demand for imaging with WGI as a multimodality system. In addition, radionuclides used for internal radiotherapy is of interest as some of these radionuclides emits single-gamma rays as well. In our future study, therefore, we will conduct WGI Compton imaging for gamma rays with energies other than 909 keV.
The system modeling could be done more accurately. The DRF derived in this paper was based on the Gaussian functions with the FWHM converted by giving perturbations to the energy information. Although the result indicated this was a sufficient approximation, the actual DRF shapes of the WGI prototype have not been measured. For validation and further improvement of the model, we will conduct the DRF measurement for developed detectors. One possible way to improve the DRF would be incorporating the Voigt function which is a convolution of a Gaussian function and a Lorentzian function with the measured parameters because the Voigt function is normally used to analyze angular resolution measure as a Compton detector performance. For the calibration of the DRF parameters, we will develop a model with a fitting function of energy and incident angle and a calibration method using measured data of a point source placed at several positions.
There is also room for improvement in data corrections. In this study, we did not apply attenuation and scatter correction methods because these effects are negligible for small subjects with high gamma ray energy. However, they are not negligible for large-sized subjects such as a human body, and imaging accuracy will be degraded. Also, these effects depend on gamma ray energy. Therefore, our future work includes imaging performance evaluation for radioisotopes with various gamma ray energies and development of attenuation and scatter correction methods for Compton imaging. We should note that a simple implementation of these correction methods is expected to have an enormous computational burden and further increase the time for reconstruction. Even at this stage, one iteration requires 2 h with two graphics processing units (NVIDIA Quadro P5000) for the mouse imaging. We should note that we have not optimized the implementation for calculation speed yet. For example, the DRF was calculated for all the voxels in each projection operation although many of the voxels do not have sensitivity. Therefore, we should develop a fast algorithm and implement the reconstruction method making corrections with an optimized projection algorithm. Also, the numbers of iterations and subsets need to be considered. We can also speed up calculations if the number of iterations can be reduced. Increasing the number of subsets improves convergence speed but image noise increases. Therefore, we are considering applying a relaxation parameter to control noise enhancement (Tanaka and Kudo 2010). These imaging parameters should be optimized by evaluating the tradeoff between noise and contrast recovery for a contrast phantom with a similar amount of radioactivity injection to an actual study.
Lastly, combined Compton imaging and PET is our final goal. For 89 Zr imaging, we reconstructed PET and Compton events separately. As both types of events originate from the same 89 Zr distribution, we can use all these events to estimate the distribution by developing an image reconstruction method combining both types. 89 Zr is one of the promising radionuclides for WGI because not only its energy characteristics but its long half-life (3.27 d) is also suitable for anti-body imaging in long-term observation. Not only for such isotopes emitting both positrons and single-gamma rays but also for positron emitters, a large number of annihilation photons can be detected as Compton events because the other photon of the pair may not be detected, and including these Compton events in image reconstruction can improve image quality (Chinn et al 2007). Therefore, combining both PET and Compton events has the potential to improve image quality by increasing sensitivity for positron emitters as well as 89 Zr. We should note that some of the Compton events can be extracted as PET events by extracting triple or quadruple coincidence events in which one side or both of the 511 keV photon pairs were scattered at the scatterer ring and absorbed at the absorber ring. We will develop a software coincidence method to extract such events and use these events for image reconstruction. Our future work also includes the development of a hybrid image reconstruction method combining PET and Compton events for WGI.

Conclusion
We developed a 3D Compton image reconstruction method, which was applied to a WGI prototype modified for small animal experiments. We measured a small rod phantom, a cylindrical phantom and a mouse with 89 Zr, which emits both a positron and a single-gamma ray of 909 keV. As a result, the small rod phantom showed Compton imaging of the prototype had spatial resolution better than 3 mm in off-center regions, and the cylindrical phantom experiment showed comparable uniformity for PET and Compton imaging of the prototype. The small animal experiment demonstrated that 89 Zr assimilated in the mouse bones was clearly depicted in both PET and Compton images. The obtained 3D mouse images are innovative in the research history of Compton imaging, and high definition Compton imaging with the image reconstruction method incorporating essential correction methods has the potential to outperform PET.