Simulation

To investigate the spectral drift related to the incident position, we constructed a model of an imaging detector based on a monolithic scintillator using the Geant4 software. The crystal used in the model is a monolithic LaBr₃(Ce) crystal, known for its exceptional properties such as density, energy resolution, light yield, and decay time. Figure 1 illustrates the structure of the detector and step-scanning process. The detector was modelled with a 51 mm×51 mm×5 mm LaBr₃(Ce) crystal coated with a Teflon foil (0.3 mm thick), which acts as the front diffuse reflective layer. The scintillator was surrounded and encapsulated by an aluminum foil (0.5 mm thick) for the specular reflective layer [17]. A glass layer (3 mm thick, reflectivity 1.5) and optical grease (0.1 mm thick, reflectivity 1.41) were placed between the bottom of the crystal and silicon photomultiplier (SiPM) array. The SiPM used in the model was the Micro-30035-TSV model manufactured by Onsemi, featuring a 16×16 array with sensitive and package areas of 3.07 mm×3.07 mm and 3.15 mm×3.15 mm, respectively [17].

Fig. 1Full-screen View

(Color online) Schematic diagram of the imaging detector simulation experiment: (a) imaging detector gridded model, (b) the response energy spectrum by irradiation at different block, and (c) detector structure

The physical properties of lanthanum bromide crystals include the emission spectrum, light yield, decay time constant, refractive index, and intrinsic resolution. The emission spectrum of LaBr₃(Ce) crystals varied with the Ce content. As the Ce concentration increased, the light yield slightly changed and the decay time decreased. In general, the decay time constant includes fast and slow time constants. As the fast decay component of the decay time constant of LaBr₃:5% Ce crystal reached 97%, the decay time constant was set to include only the fast component. The fast time constant was set to 15 ns, and the yield value was set to 1 [28, 29]. The light yield (scintillation yield) and refractive index of the crystal was set to 6.3×10⁴ MeV^-1 and 1.85 [29], respectively. In GEANT4, the emission spectrum obtained from physical experiments [30] were used to establish the emission spectrum of the crystal model. The most significant characteristic of lanthanum bromide crystals is the presence of radioactive materials ¹³⁸La and ²²⁷Ac within the crystals. The presence of these two radioactive nuclides provides lanthanum bromide crystals with inherent radioactivity. Their background energy spectrum is slightly different from the natural background spectrum. In the low-energy region (0-1800 keV), the background spectrum is mainly due to gamma transitions of ¹³⁸La (798-1436 keV). When simulating the background spectrum, we calculated the total activity of ¹³⁸La based on the volume of the detector crystal and randomly dispersed it inside the crystal (the volume activity of ¹³⁸La is 1.52 Bq·cm^-3) [31, 32], with a total activity of 19.7 Bq. Owing to the small size of the crystal, it has a lower absorption capacity for gamma rays at 789 and 1436 keV, which has minimal impact on imaging results. In the high-energy region (1800-2800 keV), the rays mainly come from ²²⁷Ac and the environmental background. As the total count of high-energy rays is small, their impact on energy spectrum measurement and imaging results is minimal.

To better simulate the transport process of optical photons, we considered the transport process of optical photons within the medium and between different media [33-35].

The effective detection area of the imaging detector composed of monolithic crystals was artificially divided into 441 blocks. Located 180 mm in front of the imaging detector, the emission position and direction of gamma rays were set using the Class PrimaryGeneratorAction to simulate a point source. The I-131 radioactive source was placed within a tungsten-shielding canister located on the incident side of the imaging detector. The tungsten- shielding canister was equipped with a collimation hole, allowing the γ-rays emitted by I-131 to be perpendicularly incident on the detector. Using a step-scanning method, each block was scanned to obtain the detector response corresponding to all positions where the rays enter the crystal. When gamma rays enter the detector from a certain block and interact with the scintillator to produce photons, these photons are converted into electrical signals by the SiPM array. The signals from 256 SiPMs were summed to serve as the detector’s energy response. Subsequently, the summed signal of the SiPM array is statistically processed to obtain the γ-ray spectrum for each incident block. Then, the peak channel of the gamma-ray spectra for each incident region was extracted. Because the event occurred near the center of the crystal, the detected light distribution was relatively complete. When gamma rays act on the edge of the crystal, the light distribution is truncated. Therefore, the peak channel of the spectrum in Block (00,00) was used as the standard peak channel.

We obtained the gamma-ray energy spectra at various detector blocks for 10 common energies. Ten distinct peaks are clearly visible in the spectra, confirming the feasibility of the established model [27]. We isolated the characteristic peak at 365 keV, as shown in Fig. 1(b). Within the detector, there is a significant spectral drift related to the incident position from Block (00,00) to Block (10,10), which intensifies with increasing energy, making the phenomenon more pronounced.

Experiment

Figure 2 shows the radiation imaging system. The system primarily consists of an encoding plate (MURA-11), imaging detector, electronic system, a Jetson Xavier NX, and display unit. The imaging detector primarily consists of a LaBr₃(Ce) crystal coupled to an SiPM array. This study utilized the SiPM array model ARRAYJ-30035-64P manufactured by Onsemi. By assembling four of these arrays, a 16×16 detector array is formed, yielding a total of 256 channels of output signals. The ARRAYJ-30035-64P consists of 3 mm×3 mm J series SiPM arranged in an 8×8 array. Four SiPM arrays were fixed to the front of the interface board through connectors on the back, leaving a 1 mm gap between each SiPM array to facilitate the packaging of the detector. The back of the interface board was composed of four pairs of connectors for the bias voltage input and initial signal readout of the SiPM arrays. All sensor substrates (cathodes) were combined to form a common I/O. After the internal packaging of the detector was completed, it was fixed to the detector housing by mounting holes at the four corners of the interface board. In the design of the interface board, owing to the lack of extensive copper plating on the board and light leakage from the connector itself, black tape was used for further light-shielding treatment of the detector. The electronic system comprising a symmetric charge division circuit, amplitude time converter circuit (ATCC), and digital signal processing circuit (DSPC) was also used. The ATCC initially employs the AD8066 operational amplifier to perform signal summation and amplification, scaling the signal to a suitable amplitude for input into a hysteresis-comparator circuit. The comparator utilized the TLV350 device, and the hysteresis comparison circuit shaped the nuclear pulse signal into a rectangular pulse signal. The 256-channel signals read out from the SiPM array detector were divided into 32-channel signals through the summing circuit and inputted into the DSPC via a flexible circuit interface for field-programmable gate array (FPGA) acquisition. The DSPC incorporates a Xilinx 7 Series core board along with a system power board. The core board captures the 32-channel rectangular pulse signals outputted by the conditioning board through its 32 I/O ports. A single-channel signal that has undergone analog-to-digital conversion was acquired by the FPGA for filtering and shaping, peak value extraction, and energy spectrum analysis. The 33-channel signal data (32-channel rectangular pulse signals and 1-channel ADC data) processed by the core board were transmitted to the Jetson Xavier NX for further processing via Ethernet. The Jetson Xavier NX is NVIDIA’s edge-computing platform, utilized for the deployment of 32-FCNN localization algorithms and maximum likelihood estimation algorithms, enabling the parallel computation of these algorithms to achieve the goal of algorithm acceleration. Ultimately, the reconstructed images and other information on the radioactive sources were displayed on the monitor. In a physical experiment analogous to simulation experiments, an I-131 source was positioned at a distance of 180 mm from the detector. A step-scanning method was utilized, in which radiation was incident perpendicularly onto the detector surface through a collimator. The detector captures the energy response at different incident positions. In the step-scanning process, a collimator with a hole diameter of 2.5 mm, overall diameter of 33 mm, and thickness of 55 mm was used to collimate the gamma rays emitted by I-131. Subsequently, the effective detection area of the detector is gridded by uniformly dividing the effective detection area of the detector into 441 blocks. The I-131 source (activity is 0.1 mCi) was placed above the collimator, marking the incident position as (x, y). Each point was measured for 4 min, and the summed signal Vxy from the SiPM array at position (x, y) was acquired through the multichannel signal acquisition board. Finally, the energy spectrum Exy at the incident position was statistically derived from the signal Vxy using the upper computer, completing the single-step scan process. Step scanning involved a total of 441 measurement points, with a total duration of 29.4 h.

Fig. 2Full-screen View

(Color online) The radiation imaging system

Correction method

In preliminary simulation studies, we extracted the peak channels of each full-energy peak in the energy range of 300-2000 keV for various blocks. Since a linear relationship exists between the channel and energy of the peak, there should also be a linear relationship between each characteristic peak in each block [27]. Therefore, in physical experiments, we can also use a functional correction matrix to improve the full-energy peak consistency. However, the functional correction matrix must be recalibrated. Thus, we investigated the linear relationship between the full-energy peak of the central Block (00,00) and other blocks through physical experiments. Because step-scanning experiments have certain requirements for the type, activity, and size of the radiation source, the laboratory can only provide I-131 that meets these requirements. Because the calibration of the functional correction matrix requires full-energy peaks of two or more energies, we selected two characteristic peaks of I-131 at 284 keV (branching ratio of 6.12%) and 365 keV (branching ratio of 81.5%) for matrix calibration. Figure 3(a) demonstrates the linear relationship between the 284 keVand 365 keV characteristic peaks of the central Block (00,00) and the peripheral Block (09,08) and Block (10,10) in physical experiments. Among them, the linear relationship formula between Block (00,00) and Block (10,10) is y_(10,10) = 0.98x − 36.29, k_(10,10)=0.98, b_(10,10)=-36.29. The linear relationships between blocks (00,00) and (09,08) were y_(09,08) = 1.02x − 38.71, k_(09,08)=1.02, and b_(09,08)=-38.71. Thus, the function correction matrix is composed of coefficients (k,b) from the linear equation. We established a linear parameter-fitting line between the spectral channel addresses of Block (i, j) and Block (00, 00) detector, with the linear relationship y_(i,j) = k_(i,j)x + b_(i,j). Here, y_(i,j) represents the channel address of the peak for Block (i, j), x indicates the channel of the peak for the central Block (00,00), k_(i,j) denotes the slope for Block (i, j), and b_(i,j) represents the y-intercept for Block (i, j). Ultimately, a function correction matrix(k_(i,j),b_(i,j)) is obtained, as shown in Eq. (1).(1)Based on the functional relationship between the incident position and energy discovered in this study, the energy-response functional correction matrix is proposed to correct the energy response and improve the consistency of full-energy peaks. First, a step-scanning method was employed to calibrate the correction matrix of the detector. This means that gamma rays were incident perpendicularly onto the detector through a collimator, and the response signals from the SiPM array of the detector were summed and statistically analyzed for energy spectra, from which the full-energy peaks and coordinate information were extracted. By using a linear function, the coefficients for the peak channel of each block relative to the peak channel of the central block were obtained, thus acquiring the complete energy-response functional correction matrix of the detector. Ultimately, the coordinate information of the interaction between the rays and the detector is used. Based on the principle of the “look-up table’’ method, the correction coefficients corresponding to the coordinates in the energy-response functional correction matrix were found. Therefore, the peak channel of each interaction event at different locations and the energy response of the entire detector were corrected, as shown in Fig. 3(c). As seen in the figure, E_(i,j) is the uncorrected peak address, (k, b) is the functional correction matrix, and is the corrected peak address.

Fig. 3Full-screen View

(Color online) (a) Relationship between peak positions at different locations and the central position. (b) Neural-network architecture, and (c) correction method

This imaging system utilizes a ray localization algorithm based on a neural network to obtain positional information on the rays [35]. The principle is that when rays interact with the detector, the SiPM at the interaction point receives the strongest pulse signal intensity. As the scintillation light propagates within the crystal, the signal strength of the SiPM at other locations gradually decreases. Therefore, the position of the rays can be determined based on the distribution of light intensity across the SiPM array. In this imaging system, the 256-channel SiPM signals were summed and converted into 32-channel signals, which were used as input data for a fully connected neural network to perform neural-network ray localization to obtain positional information. The structure of the neural network model includes one input layer, two hidden layers, and one output layer. The input layer consists of 32-channel signals. The output layer consists of the coordinates (X, Y) of the interaction position between the incident ray and crystal. The rectified linear unit activation function was used. On this basis, the intrinsic spatial resolution of the radiation imaging system can reach 2.67 mm. The modified uniformly redundant array (MURA) can be represented by a matrix Bij, where r denotes the number of rows (or columns) in the matrix, which is a prime number. Bij=1 and Bij=0 indicate an open and closed aperture, respectively [36].(2)where (3)The structure of the nested MURA-coded aperture was derived from the MURA (N × N) by expanding the MURA-coded aperture in a centrally symmetric manner to nine times its original size, resulting in the loop MURA (3N × 3N). Subsequently, a coded aperture of size (2N-1)×(2N-1) at the central position was selected to serve as the nested MURA [36]. After the interaction of radiation with the crystal, a neural-network localization algorithm was employed to acquire the encoded image on the detector.

Correction results

To better compare the uncorrected and corrected results from both the simulation and physical experiments, the measurement object in this study for both simulation and physical experiments is the characteristic gamma rays emitted by I-131. The correction results will be discussed in terms of the consistency of uncorrected and corrected peak positions, energy resolution, encoded images, and reconstruction image quality.

Given the limited variety of radiation sources in the laboratory, to more rigorously compare the simulation and experimental results, the spectral channel corresponding to the 365 keV peak was isolated from the simulation dataset. This extraction provided the 365 keV peak channel across various blocks for subsequent analysis, as shown in Fig. 4(a). The correlation between the 365 keV peak channel and its spatial positions indicated that the displacement is minimal at the center and increased toward the periphery, suggesting a positive correlation between the displacement magnitude and proximity to the central axis. The peak displacement is minimal at the center and gradually increased toward the periphery. This phenomenon is caused by the interaction of gamma rays with the crystal lattice, leading to the ionization and excitation of lattice atoms, followed by photon emission during de-excitation. The interaction of photons with the crystal includes the processes of absorption, refraction, and reflection, with backscattering playing a more significant role in peripheral areas. Photons traversing the crystal are modulated by the crystal structure and physical properties, leading to photon flux attenuation. This, in turn, results in the spatial heterogeneity of the detector’s intrinsic response characteristics, leading to position-dependent energy responses to gamma rays.

Fig. 4Full-screen View

(Color online) (a) the uncorrected distribution map of the I-131 peak in the simulation(@365 keV), (b) the uncorrected distribution map of the I-131 peak in the experiment(@365 keV), (c) the corrected distribution map of the I-131 peak in the simulation(@365 keV), (d) the corrected distribution map of the I-131 peak in the experiment (@365keV), (e) uncorrected and corrected I-131 energy spectra in the simulation, and (f) uncorrected and corrected I-131 energy spectra in the experiment

Figure 4(b) illustrates the distribution of the 365 keV characteristic peaks in the experiment. The experimental data closely align with the results obtained from MC simulations. The variation in pulse-height channel addresses were pronounced owing to the different incident positions of the radiation. During the interaction of gamma rays with the crystal, ionization and excitation of crystal atoms occurred, and the de-excitation process emitted photons. These photons may be absorbed, refracted, and reflected by the crystal, with the edge positions significantly affected by backscattering. The propagation of photons within the crystal is influenced by its structure and physical properties, causing light loss. This results in inconsistent energy responses of monolithic crystal detectors. The overall peak position indicates that the closer to the center position, the smaller the spectral drift related to incident position. The closer to the peripheral position, the more severe the spectral drift related to incident position, which is similar to the simulation results. In the experiment, the degree of the peak position offset was more severe, reaching 84 channels. Although SiPMs exhibited a certain degree of temperature drift, the entire experimental process was conducted at room temperature, with the signal processing circuit board and FPGA chip maintained at a constant temperature. Therefore, the spectral drift caused by temperature can be considered negligible owing to the difference in the uniformity of the crystal itself. In the simulation, the LaBr₃(Ce) crystal was isotropic, and properties were exactly the same everywhere within the crystal. The LaBr₃(Ce) crystals used were not uniform, which was manifested in the specific growth direction of the crystals and differences in Ce doping concentrations at various locations within the crystals. These differences will cause the inhomogeneity of the anisotropic transmission of photons within the crystal and the difference in the number of photons emitted when depositing scintillation events of the same energy. These inconsistencies will affect the amplitude of the pulse signal of the detector, thereby influencing the energy resolution. The differences in crystal packaging can also lead to this result. In the simulation, not only is the crystal itself absolutely uniform, but the packaging of the crystal surface was consistent everywhere. However, it is difficult to achieve this consistency experimentally. The differences brought about by the packaging of the LaBr₃(Ce) surface will prevent the light transmission process from being completely uniform and symmetrical. When charged particles are incident from different positions on the crystal surface, geometric factors affect the transmission of photons within the crystal, eventually resulting in different peak positions.

From the consistency of full-energy peaks, the energy-response functional correction matrix was obtained separately from the simulation and physical calibration experiments. The energy response for both simulated and physical experiments were corrected, resulting in a significant improvement in full-energy peak consistency. Figure 4(c) illustrates the corrected 365 keV peaks for various incident blocks in the simulation. Figure 4(d) shows the same for the experimental results. It can be observed that the peaks in the edge blocks were significantly corrected in the simulation and experiment. The consistency of the full-energy peak was significantly enhanced, essentially converging to uniformity. In addition, within the energy range of 365-1332 keV, the results of the simulation are well-fitted to the energy-response functional correction matrix [27]. The energy-response functional correction matrix was designed to channel the 284 and 365 keV energies of I-131 in the experiment. Consequently, the peak channel error was corrected to zero channels in the physical experiment.

Owing to inconsistent responses to uncorrected peaks, there is a variation in the energy resolution of the detector. When incident gamma rays are predominantly absorbed by the edge blocks of the detector, the energy resolution of the detector decreased. Conversely, when the central block of the detector primarily absorbs incoming rays, the energy resolution improved. Thus, enhancing the consistency of the detector’s full-energy peak can improve the energy resolution to a certain extent. The energy-response functional correction matrix was obtained separately from the simulation and experiment to correct the energy response, resulting in improved energy resolution. The energy resolution at the 365 keV energy, uncorrected and corrected in the simulation and experiment, are shown in Figs. 4(e) and (f), respectively. The full-width at half maximum (FWHM) of the corrected gamma-ray spectrum was narrowed, the peak count was obviously increased, and the energy resolution was improved. In the simulation, the FWHM of the characteristic peak at 365 keV was reduced from 33 to 28 channels, increasing the energy resolution by 15.2%. In the experiment, the FWHM of the characteristic peak at 365 keV was reduced from 53 to 38 channels, increasing the energy resolution by 28.3%. The improvement of the energy resolution of the imaging detector plays a key role in selecting the energy window data and improving the quality of the reconstructed image.

Image reconstruction

In imaging applications of coded-aperture gamma-imaging systems, data is selected based on the energy window. Owing to the inconsistent response to uncorrected peaks, the detection efficiency and the loss of real events will be reduced, which will affect the image quality. In more severe cases, this can lead to incomplete projected images and ultimately failing to form an image, particularly when the incident rays are absorbed by the edge block of the detector. Therefore, it is necessary to correct the full-energy peak consistency in the effective detection area of the detector to improve the imaging quality. In this study, both the simulation and physical experiments utilized a radioactive source of I-131(d_1/2=8.02 d) with an activity of 0.1 mCi. The center of the I-131 point source, encoding board collimator, and imaging detector were aligned horizontally along the same axis. The distance between the I-131 point source and detector was 80 cm, and the distance between the encoding board collimator and detector was 5 cm.

Figures 5(b) and 5(c) represent the uncorrected and corrected encoded images from the experiment, respectively. In Fig. 5(b), data loss was observed in the upper-right and lower-left corners, with a more significant loss in the lower-right corner. This was attributed to the loss of true events in the upper-right and lower-left blocks, as shown in Fig. 4(b), leading to imprecise energy window data. The uncorrected encoded images in the physical experiment were severely hampered by data loss, rendering image reconstruction infeasible. In addition to the event positions causing peak inconsistencies, the manufacturing process of SiPMs cannot guarantee identical performance for each SiPM. The response of the SiPM array, signal acquisition, and readout were not perfectly synchronized, which further aggravated the inconsistency of full-energy peaks and reduced the amount of effective encoded image data. The corrected encoded images are more complete, enabling a clearer reconstruction of the images. When the radioactive source is located at the center of the useful field-of-view (UFOV) of the image detector, the most valid information is collected in the central area of the detector. An incomplete encoded image will not lead to a failure in image reconstruction, and it is also possible to accurately reconstruct the position of the radioactive source. In this case, inconsistencies in full-energy peaks can introduce artifacts into the reconstructed image, reducing the signal-to-noise ratio (SNR) and affecting the quality of the reconstructed image. However, when the radioactive source is located at the edge of the UFOV of the image detector, most of the valid information is concentrated in the peripheral area of the detector. The existence of inconsistent full-energy peaks can cause severe loss of valid information, leading to deviations in the position of the radioactive source in the reconstructed image, or even failure to image. Figures 5(d) and (e) represent the uncorrected and corrected reconstructed images from the experiment, respectively. The SNR of the uncorrected reconstructed image is only 2.38. The corrected peak channel eliminated data that does not belong to the energy window and improved the efficiency of the detector. Simultaneously, it improves the energy resolution of the imaging detector. The SNR is significantly enhanced, from 2.38 to 5.37, an increase of 125.6%. It can be seen from Fig. 5 that the energy resolution of the uncorrected imaging detector is poor, resulting in the data not belonging to the energy window recognized as valid data. This results in a low SNR and the presence of artifacts in the reconstructed image. The SNR of the corrected reconstructed image is improved, reducing or eliminating artifacts, and enhancing the quality of the reconstructed image.

Fig. 5Full-screen View

(Color online) (a) Location coordinates of I-131 source. (b) Uncorrected encoded images. (c) Corrected encoded images. (d) Uncorrected reconstructed images, and (e) corrected reconstructed images