A. Principle of Compton Imaging

As an electronic collimation imaging method without collimator blocking particles, Compton imaging can provide a better efficiency than traditional mechanical collimation techniques as it can specify the incident direction of the gamma-ray to the surface of a cone using multiple-interaction events [7]. The kinematic reconstruction of Compton scattering is illustrated in Fig.1. As an example, consider a two-sites full energy deposition event; the energy of the incoming photon, $E_0$, and the scattering angle, $\text{cosθ}$, can be determined by:

Fig.1.Full-screen View

(Color online) Example of a Compton scattering interaction and Compton back-projection cone.

where $E_1$ is the energy of the recoil electron, $E_2$ is the energy of the scattered photon, $m_{\text{e}}$ is the rest mass of the electron, and c is the speed of light [8].

With the scattering angle determined, all possible directions of the incoming photon form a surface of a cone; the axis of the cone is given by the scattered photon direction and the opening angle is equal to the scattering angle, as shown in Fig.1.

B. Principle of 3D position sensitive CZT detector

In order to realize Compton imaging, we need to record the deposited energies and the positions of multiple-interaction events occurring in the active volume. The 3D CZT detector can fully meet these demands. It has pixel electrodes on the anode side and a conventional planar electrode on the cathode side. Similar to a time projection chamber (TPC) [9], the (x, y) coordinates of interactions are obtained from the location of the pixel anodes where electrons are collected, and the interaction depths (z-coordinate) between the cathode and the anode are obtained from electron drift time [5].

In our study, we focus on the two-site interactions full energy deposition events consisting of Compton scattering followed by photoelectric absorption. As illustrated in Fig.2, an incoming photon undergoes a Compton scatter at position $\overrightarrow{r_1}$, creating a recoil electron of energy $E_1$ and then the scattered photon undergoes a photoelectric absorption at position $\overrightarrow{r_2}$ and depositing energy $E_2$. By measuring the positions and energies of the two interaction sites in the 3D CZT detector, the Compton cone can be reconstructed by the Compton formula (Eq.2). As increasingly more events are recorded, many back-projection cones are reconstructed and the intersection of these cones can be used to verify the localization of the radioactive sources.

Fig.2.Full-screen View

(Color online) Principle of Compton imaging by a 3D position sensitive CZT detector.

A. Sequence reconstruction

In the actual Compton imaging process, we usually do not have a priori knowledge regarding the radioactive sources in the environment. An energy window will be set according to the measured energy spectrum, and all events with energy within this window will be used to reconstruct the source image. These are termed full energy deposition events implying that they end with photoelectric absorption.

For those events with energy that falls out of the energy window, termed a partial energy deposition event, the incident gamma-ray energy can only be assumed to be the sum of several individual deposited energies (Eq.1). In this case, the usage of partial energy deposition events will also cause a wrong back-projection cone (Fig.3). Therefore, only the back-projection cone from full energy deposition events with the right sequence can pass through the actual source position.

For simplicity, we only considered the sequence reconstruction method for two-site full energy deposition events. As the Compton scattering angle increases from 0 to 180 degrees, the energy of the recoil electron increases gradually. The maximum deposited energy in the first scattering can be calculated by Eq.3, which is known as the Compton edge [11] in a typical gamma-ray energy spectrum; a gamma-ray photon cannot lose more energy than this energy edge. In other words, if any one of the two energies detected in the detector is greater than the Compton edge, that energy must belong to the photoelectric absorption interaction.

If the energies deposited in both interactions are smaller than the Compton edge, two techniques can be applied to determine the right sequence of the interactions; termed the "simple comparison" and "deterministic algorithm" methods [12]. The first method is based on the knowledge that the higher one of two energies deposited in a detector is more likely to belong to the scatter interaction for gamma-rays with energies greater than 400 keV. The second method is based on the Klein-Nishina formula [13]; the probability of Compton scattering followed by photoelectric absorption will be computed for two possible sequences, and the sequence with the higher probability is considered to be the correct one. The latter method can also be used for the sequence determination of events with more than two interactions.

B. Far-field approximation

During the reconstruction process, the vertices of each back-projection cone should be the point where scattering occurs and these vertices are randomly distributed in the detector volume. However, for most application scenarios, the detector size is quite small compared with the distance between the radioactive sources and detectors. Therefore, the spatial differences between these vertices can be ignored and all reconstruction cone vertices are approximated at the center of the detector. Then the intersection of these cones will only convey the directional information of the sources because all cones come from the same vertex [Fig.5(a)].

Fig.5.Full-screen View

(a) Far-field imaging only conveys directional information. (b) Near-field imaging can localize sources in 3D space.

A two-dimensional coordinate system is now built for the image space, as shown in Fig.4. By assuming that there is an image sphere surface at infinity, the detector can be considered as a point at the origin of the sphere. All back-projection cones intersect the sphere into a ring. Then the complicated process of finding the intersection of Compton cones in a 3D space is simplified by finding the intersection of Compton rings in a 2D sphere; this process is known as the "far-field approximation". As a result, the calculation process is greatly simplified but the radial information is lost [Fig.5 (a)], and every direction on the reconstruction cone can be calculated easily through the angular formula (Eq.4):

Fig.4.Full-screen View

In a far-field approximation, the calculation is greatly simplified.

where θ_c is the Compton angle calculated from the deposited energy; θ₂ and φ₂ are the azimuthal and polar angles of the back-projection cone axis that can be calculated from the positions of the scattering interaction and absorption interaction; and θ₁ and φ₁ represent the possible direction of source.

However, if the source is in the vicinity of the detector, the spatial differences between these vertices of different back-projection cones must be taken into account [Fig.5 (b)]. The separation distance of each of the individual Compton cone vertices provide enough parallax for the localization of a radioactive source in 3D space. This is termed "near-field imaging" [14], however, for a single stationary detector, most scenes are "far-field approximations".

C. Simple back-projection

The simple back-projection method is the most straightforward algorithm. During the reconstruction process, each event recorded in the detector is reconstructed in its own individual image by giving the same weight to all possible source directions that are then added to the overall image. Thus, the reconstruction is in real-time and the overall image can be updated for each event. Finally, the true source direction will be enhanced if we have a fairly good measurement.Fig. 6

Fig.6.Full-screen View

Simple back-projection rings are added to the overall image.

D. Filtered back-projection

Even with the perfect energy and position resolutions, correct sequence, and without Doppler broadening effects [7], there is still a significant blur in the SBP result because the entire circle has a probability distribution for the event originating from a certain direction. Parra proposed a FBP algorithm in the spherical harmonics domain [15]; an ideal point source will generate a simple back-projection result, as in Eq.5, and is clearly not a theoretical delta function.

Here, $\omega$ is the angle between the image pixel and source directions, $\left(\stackrel{\rightharpoonup }{\Omega }-{\stackrel{\rightharpoonup }{\Omega }}^{\text{'}}\right)$, and $K\left(\theta \right)$ is the probability density function of the photon scattered in the direction $\theta$ defined by the Klein-Nishina formula. Equation (5) is the point spread function (PSF) of the system. If the detector geometry has spherical symmetry, the PSF will be shift-invariant, which means that the PSF remains invariant to the photons incident from all directions.

Suppose the SBP image is $g^{\prime }\left({\stackrel{\rightharpoonup }{\Omega }}^{\text{'}}\right)$, the true source image is $g\left(\stackrel{\rightharpoonup }{\Omega }\right)$, and their relationship must obey Eq.6, in which $\text{*}$ is the convolution operation. According to the Fourier convolution theorem, the Fourier transform of a convolution is the pointwise product of the Fourier transform, leading to Eq.7. The Fourier transform for an arbitrary function, $g\left(\stackrel{\rightharpoonup }{\Omega }\right)$, defined on the 4 $\text{π}$ sphere can be transformed into the spherical harmonics by Eq.8 [15].

Here, $Y_l^m\left(\stackrel{\rightharpoonup }{\Omega }\right)$ are the spherical harmonics on the 4 $\text{π}$ sphere, $Y_l^{m*}\left(\stackrel{\rightharpoonup }{\Omega }\right)$ are its complex conjugates, and $P_l^m\left(\text{cos}\theta \right)$ are the associated Legendre polynomials.

We can then obtain the Fourier transforms of the SBP image and shift-invariant PSF, and the FBP image can be calculated according to Eqs.7 and 9.

A Compton camera with a 121 pixels CZT detector was simulated by Geant4 for the numerical experiments. The size of the detector was 20 mm by 20 mm by 15 mm, the pixel pitch was 1.72 mm, and a steering grid between the pixels was used in the simulation. Two separate point sources were placed 50 cm from the detector center. The result of the SBP [Fig.7 (a)] and FBP [Fig.7 (b)] images of the simulated two-point source are shown below.

Fig.7.Full-screen View

(Color online) (a) SBP and (b) FBP images of a simulated two-point source.

E. Maximum likelihood expectation maximization

The maximum likelihood expectation maximization (MLEM) is an iterative algorithm that determines an estimate image intensity distribution that maximizes the probability of observing the actual detector count data over all possible density distributions [16]. Although the MLEM is not intended for the reconstruction result of a point source for a specific detector system, our goal is to compare it with the results obtained with FBP. The iteration is performed by the following equation:

where ${\lambda }_j^n$ is the estimated density distribution at step n, $g_i$ is the observed measurement, $T_{ij}$ is the system matrix describing the probability that a photon from direction j is observed as a measurement i, $s_j$ is the sensitive matrix describing the probability that a photon from direction j is recorded by the detector. It can be calculated using the Monte Carlo method proposed by Xu [7]. This algorithm has several characteristics: the likelihood strictly increases at each step; non-negative results will be obtained if the iterations are started from a non-negative initial value; and the total number of estimated counts is always equal to the total number of observed counts at each step. The iterations can be stopped by the accuracy requirements or computation time limits.

In the MLEM algorithm iteration process, it is necessary to sum the probability over all possible measurement outputs, i, but there are too many possible detector responses. It is unrealistic to pre-calculate and sum the whole system matrix. An alternative method is to compute only the probability of those measured events and this data storage mode is called the list mode [17].

The system matrix $T_{ij}$ can be modeled by the attenuation probability according to the path that a gamma-ray travels in the detector material (Eq.12). Fig.8 shows the process of a gamma-ray from direction $\left[{\theta }_1,{\phi }_1\right]$ incident from A, scattering at point B, and finally absorbed at point C. The corresponding reaction section in the formula can be found in the NIST database [18].

Fig.8.Full-screen View

Example of the path for a two-interaction event in the detector.

A. CZT detector and ASIC

A Redlen CZT crystal was fabricated by Nuctech Company Limited. The detector is 10 by 10 by 10 ${\text{mm}}^3$ and has 4 by 4 pixelated anodes with a pitch of 2.3 mm and a planar cathode. The structure of the detector can be accessed at Ref.11. We used the VASTAT ASIC [19] designed by IDEAS to read out the amplitude and timing information of the CZT detector. Each channel has a charge sensitive pre-amplifier and shaper circuit for the amplitude measurement. Besides the slow shaper, each channel also has a fast shaper, followed by a discriminator and time-to-amplitude converter (TAC) for measuring the electron drift time.

B. DAQ and software

The ASIC module is driven by a data acquisition (DAQ) system. The DAQ board is equipped with an FPGA that can generate and send the readout clock signals to the ASIC and also convert the output of the ASIC voltage signals to a digital signal by the A/D converters. The DAQ board is also used as the controller interface between the DAQ software and ASIC module. A personal computer (PC) is used to communicate with the DAQ board through Ethernet. The PC serves as a master control device and storage media. A software application based on the LabVIEW has been developed to control the acquisition process and allows the user to configure the ASIC.

Fig.10 shows the main windows of the DAQ software. The system setting window can set the thresholds of the anode and cathode, and configure parameters of the ASIC. The real-time display window monitors each received event and provides a brief on-line analysis for a quality check of the collected data. The received data can be saved in a binary-form file for offline analysis. The offline data analysis program was written in C++ and based on the ROOT [20] with the function of event reconstruction and Compton imaging.

Fig.10.Full-screen View

Main windows of the DAQ software.

A. Calibration and spectra performance

For multiple-pixel events, the gamma ray energy is estimated by summing multiple-pixel signals together. Because of the material non-uniformity in the detector, and the baseline offset and gain variation in the ASIC channels, calibration should be carried out to obtain accurate energy.

The detector was biased to -700 V at room temperature. We calibrated the detector with the 276, 306, and 356 keV lines from Ba-133 and the 662 keV line from Cs-137. The photopeak positions (ADC channels) as a function of gamma ray energies are shown in Fig.11. A linear function is fit to the curve. The gain is determined from the slope of the fitting line and the intercept indicates the baseline offset (Fig.12). All 16 pixels show good linearity and have almost the same gain. However, the baseline offsets vary significantly from channel to channel.

Fig.11.Full-screen View

Relationship between the ADC channels and energy for pixel #3.

Fig.12.Full-screen View

Gain and offset of 17 channels; channel #0 is the cathode channel and the remaining 16 are anode channels.

Depth sensing and correction techniques were used to reconstruct the positions of the interaction and deposited energies of each gamma photon. The depth of the interaction was obtained by the C/A ratio or electron drift time. The whole detector was divided into $16\times 15$ voxels after we used 15 bins along the depth dimension. The photopeak position in each voxel-based Cs-137 spectrum of single-pixel events was used to obtain a gain factor matrix, then the energy was corrected according to this gain factor matrix for each event. A relationship between the electron drift time and C/A ratio was established to correct the interaction depth for multiple-pixel events. After correction, the detector shows good spectra performance. The single-pixel events spectra of 16 pixels are shown in Fig.13. The energy resolution (FWHM) of most pixels for the 662 keV line is better than 2.0%. An electronics noise of 8.2 keV was measured with a test pulse, which is the main limitation for the detector energy resolution. Efforts will be made to reduce the electronic noise in the future; a better power isolation will be performed and the power decoupled to make it clean, with a separation of the digital and analog parts of the circuit to reduce crosstalk.

Fig.13.Full-screen View

Single-pixel events spectra for Cs-137 gamma-ray source.

B. Imaging reconstruction results

A 10 μCi Cs-137 point source was placed above the cathode surface of the detector at a distance of 11 cm. The trigger thresholds for the anode and cathode were approximately 150 keV and 50 keV, respectively. Only these signals with amplitudes that exceed the threshold will trigger the ASIC to record this event. Forty thousand events were acquired. An offline Compton reconstruction was performed to image the source, and is described as follows.

Firstly, we selected the data of two-pixel events with total deposited energies within the energy window 620–700 keV, leaving only 1257 events. In addition, events with energies in two adjacent pixel events were removed because charge sharing events cannot be easily distinguished from Compton events, excluding a further 1001 events. After the selection, there were only 265 events out of total 40 thousand events that could be used to reconstruct Compton cone.

Secondly, a sequence reconstruction was performed. The interaction with the higher energy deposition was selected as the first interaction point (i.e., scatter point) and the scattering angle was calculated.

Thirdly, as the source-to-detector distance of 11 cm is much larger than the detector size of 1 cm, the far-field approximation was made to simplify the calculation. The simple back projection (SBP) method was used to project the possible incident direction (i.e. Compton cone) back into the imaging space. Intersections of the cones yielded the direction of the source. The obtained reconstructed images from a sample of four events and all 265 events are given in Figs.14 (a) and 14(b). The MLEM method was also performed and a better angle resolution was achieved after 20 iterations [Fig.14 (c)]. The number of iterations was chosen by stopping before the image become too noisy.

Fig.14.Full-screen View

(Color online) Reconstruction images. (a) SBP result with four events, (b) SBP result with 265 events, and (c) MLEM result with 265 events after 20 iterations.

Fig.15 shows a cross-section along the phi-axis of the images shown in Figs.14 (b) and 14(c) through the center of the source. The distribution obtained with the back-projection algorithm shows large tails, which are not present in the distribution obtained with the MLEM algorithm. The angular resolutions of the point source in Fig.14 (b) and (c) are approximately 36 and 8.3 degrees FWHM respectively.

Fig.15.Full-screen View

Cross-section along the image plane for the azimuthal direction.

Thereafter, a two cesium source measurement was performed under the same system conditions, such as the trigger threshold and noise level. Two sources with identical intensities were placed in the vicinity of the detector with polar and azimuthal angles of [60, 90] and [60, 180], respectively, and the far-field approximation condition was satisfied. Based on the 986 imageable events of the total 50000 events acquired during the measurement, both the SBP [Fig.16 (a)] and MLEM [Fig.16 (b)] reconstruction methods were employed, and the angle resolution was approximately 39 and 9.5 degrees FWHM, respectively. Fig.17 shows a cross-section along the phi-axis of the images shown in Figs.16 (a) and 16(b) through the center of the source.

Fig.16.Full-screen View

(Color online) Reconstruction images. (a) SBP and (b) MLEM with 968 events.

Fig.17.Full-screen View

(Color online) Cross-section along the image plane for the azimuthal direction.