Experimental setup

Figure 1(a) illustrates the experimental setup. An ²⁴¹Am-Be radioactive neutron source with an activity of 0.3 Ci was placed at the center of a paraffin box and neutrons were led through a channel with a diameter of 10 cm. As shown in the inset in Fig. 1(a), an object containing Cd and NaCl was placed along the path of the neutron beam. Four 1 mm thick Cd plates were positioned at the four corners whereas five 5 mm thick NaCl plates were placed in the remaining regions. When the target nuclides (¹¹³Cd and ³⁵Cl) captured thermal neutrons through ³⁵Cl (n, γ) ³⁶Cl and ¹¹³Cd (n, γ) ¹¹⁴Cd reactions and entered the excited state, various gamma rays were produced from their de-excitation. The inset shows an image of the sample and illustration of thermal neutron capture reactions. Figure 1(b) and (c) show the partial decay schemes of ³⁶Cl and ¹¹⁴Cd, respectively, with high thermal neutron-capture cross-sections. The gamma rays used for image reconstruction were 558 keV (Cd) and 1164 keV Cl) because of their large cross-sections. Various collimators with random patterns of blocked pixels were placed behind each sample. The characteristic gamma rays were spatially modulated after passing through the collimator mask and collected using a hyperpure germanium (HPGe) detector.

Fig. 1Full-screen View

(Color online) (a) Schematic diagram of the neutron-induced gamma ray imaging setup. (b), (c) Decay scheme of ³⁶Cl and ¹¹⁴Cd for partial gamma rays, respectively

Single-pixel imaging method

In neutron-induced gamma-ray imaging, the net peak area of prompt gamma rays emitted from a sample can be calculated using the following equation [11]: $A=\frac{m}{M}\Phi N_{\text{A}}\sigma \epsilon t,$ (1) where A is the net peak area, m (g) and M (g ⋅ mol^-1) are the mass and molar mass of the nuclide, respectively, Φ (cm^-2·s^-1) is the average thermal neutron flux within the sample, NA is Avogadro’s number, σ (cm²) is the partial thermal neutron capture cross-section, ε is the detection efficiency, and t (s) is the measurement time.

The basis of this technique is to use various masks and detect the signals of correlations between the masks and the sample. An image can be reconstructed by multiplying each mask by the corresponding signal. For an image with the total number of pixels N = P₁× P₂, Eq. (1) can be expressed as $A={\displaystyle \sum _{n=1}^N}{\Phi }_n{\epsilon }_n{m}_n\frac{N_{\text{A}}\sigma }{M}t,$ (2) where mn is the mass of the nuclide in the nth pixel, Φn is the thermal neutron flux in nth pixel, and εn is the detection efficiency for nth pixel. For each pixel, it becomes an isotropic volume gamma ray source with pixel dimensions after thermal neutron irradiation. These gamma rays move through and interact with samples. Because of these interactions, some gamma rays lose energy and disappear, whereas others exit the sample and traverse the collimator to the HPGe detector, where they are counted. In this study, the ratio of the counted gamma rays to the number of gamma rays emitted from the pixel volume source is defined as the pixel detection efficiency.

Inspection of Eq. (2) reveals that the parameters M, N_A, Φn, σ and t are constants for a specific sample and the given measurement conditions. The measured net peak area depends on mn and εn. If the aim is to use only L (L<N) measurements for image reconstruction, the principle of the proposed method can be summarized as follows: $A=\left[{\epsilon }_n\right]\left[m_n\right],$ (3) where [A] is a L×1 column vector, [mn] is the image with N pixels ordered in a N×1 vector, [εn] is a L×N detection efficiency matrix, whose l_th (1≤l≤L) row corresponds to the l_th measurements. The number of pixels in an image depends on the efficiency matrix.

In this study, the maximum likelihood expectation maximization (MLEM) algorithm is used to reconstruct the image. The MLEM algorithm, which is valid for data with a Poisson distribution, has been widely used for image reconstruction in computed tomography (CT), positron emission tomography (PET), and single-photon emission computed tomography (SPECT) [33-35]. The gamma rays induced by neutrons are Poisson distributed, and compared with traditional numerical algorithms, MLEM can achieve high-resolution image reconstruction through iterative calculations, especially for cases with low SNR.

The MLEM algorithm is used to calculate the most likely distribution [m] from the measured [A] by assuming that the values Ai of [A] are Poisson distributions. This algorithm is expressed as follows [36]: $m_n^{k+1}=\frac{m_n^k}{{\displaystyle {\sum }_l{\epsilon }_{ln}}}{\displaystyle \sum _i{\epsilon }_{ln}}\frac{A_i}{{\displaystyle {\sum }_n{\epsilon }_{ln}}{m}_n^k},$ (4) where $m_n^k$ and $m_n^{k+1}$ are the reconstructed values for the pixel n after k and k+1 iterations, respectively. εln is the detection efficiency of pixel n in l_th measurement.

Coded collimator

For neutron-induced gamma ray imaging technology based on a small-yield neutron source, an object can be considered as a gamma ray source with low activity after irradiation. To obtain a satisfactory imaging performance, the ratio of the open squares was set to 50%. The random coding method was used to design the collimator patterns. A set of 30 collimators made of iron were used in this study. The schematic of a typical collimator pattern is shown in Fig. 2(a). The iron pixels shield the gamma rays, whereas the opening pixels are transparent to gamma rays. The collimator dimensions were 6 cm×6 cm×10 cm ($X\times Y\times Z$). For each pattern, half of the area was blocked and the size of each square opening was set to 1 cm×1 cm.

Fig. 2Full-screen View

(Color online) (a) Schematic diagram of one of 30 collimators. (b) Monte Carlo simulations used to calculate the detection efficiency matrix [εn]. (c) Patterns of 30 collimators

In order to solve Eq.(4), the prior detection efficiency matrix [εn] is obtained. The efficiency of a certain energy depends on the HPGe detector performance and the measurement conditions, including the distance between the sample and detector, the geometry, and the materials of the sample. For a small sample, the efficiency curve can be experimentally calibrated using reference gamma ray sources. When the sample is large and cannot be viewed as a point source, the gamma ray self-absorption effect must be corrected. In this study, a Monte Carlo simulation was performed to obtain the efficiency matrix because it is difficult to determine the value experimentally. An image containing 30×30 pixels was reconstructed from 30 measurement values to obtain an image with pixel dimensions of 2 mm×2 mm. The efficiency matrix is presented in Fig. 2(b). Each row of [εn] is associated with a collimator pattern. When the Cd and Cl in each pixel capture thermal neutrons, they can be considered gamma ray sources with energies of 578 and 1164 keV, respectively. Thus, for each pixel, a volume gamma ray source with isotropic and mono-energy pixel dimensions was defined and placed at the corresponding position. The energy deposition in the HPGe detector was calculated to obtain its detection efficiency. The patterns of the 30 collimators are shown in Fig. 2(c). The efficiency was associated with the collimator pattern, as the iron pixels shielded the gamma rays whereas the opening pixels were transparent to gamma rays.

Benchmark of neutron source

The thermal neutron flux at the sample position was evaluated using an In foil activation method. Foils covered with Al and Cd were irradiated at the sample position. After cooling, the foil was measured using a LaBr₃ detector in the lead chamber. The reaction rate R can be expressed as follows [37]: $R=\frac{AM\lambda }{N_{\text{A}}mfp\epsilon \left(1-e^{-\lambda t_1}\right)e^{-\lambda t_2}\left(1-e^{-\lambda t_3}\right)},$ (5) where A is the net count of delay γ ray peak, M (g·mol^-1) is the atomic mass of isotope, λ (s^-1) is the decay constant, N_A is Avogadro’s number, m (g) is the mass of the foil, f is the isotopic abundance, p is the emission probability, ε is the detection efficiency, t₁, t₂ and t₃ (s) are the irradiation time, cooling time and measurement time, respectively. Then, the thermal neutron flux ϕ_th and epithermal neutron flux ϕ_epi can be obtained as follows [37]: ${\varphi }_{\text{th}}=\frac{1}{G_{\text{th}}g{\sigma }_{\text{In}}}\left[R_{\text{Al}}-R_{\text{Cd}}\left(1+\frac{g\sigma }{G_{\text{res}}I}{f}_1+\frac{\sigma w^{\prime }}{G_{\text{res}}I}\right)\right],$ (6) ${\varphi }_{\text{epi}}=\frac{R_{\text{Cd}}}{I},$ (7) where G_th and G_res are the thermal and epithermal self-shielding factors, respectively. ${\sigma }_{\text{In}}$ is the thermal neutron reaction cross-section. In the thermal energy range, the cross section varies inversely with the neutron speed (1/ν). g is the correction factor that accounts for departures from the ideal 1/ν cross-section. I is the resonance integral cross section, g is the correction factor caused by the neutron energy deviation of 1/ν, which is related to the neutron temperature, f₁ is the correction caused by the energy from 5 kT (the lowest energy of epithermal neutrons is usually taken to be equal to 5 kT) to the epithermal neutrons, w’ is the correction caused by the energy from 5 kT to E_Cd. The detailed values are listed in Table 1.

A delayed gamma ray 1293 keV of ^116mIn was used to calculate the thermal neutron flux. The measurement conditions of two foils and some parameters are listed in Table 1. The thermal and epithermal neutron fluxes are 40.2 ± 1.3 cm^-2s^-1 and 0.52 ± 0.01 cm^-2s^-1, respectively. The uncertainties are mainly derived from the counting statistics of the gamma rays. Compared with large facilities, the thermal neutron flux is very low. The MCNP5 software combined with the ENDF/B-VI cross-section library was used to calculate the thermal neutron flux. The MCNP model of the neutron source was developed in our previous study [38]. The neutron yield of the ²⁴¹Am-Be source was also evaluated. The thermal neutron flux (<0.5 eV) in the In foil (ϕ_MCNP) was calculated by using F4 tally, and the value was 8.29×10^-4 cm². Then, the neutron yield can be calculated using ${\varphi }_{\text{th}}/{\varphi }_{\text{MCNP}}$. The result was 4.82×10⁵ s^-1. In addition, the thermal neutron flux distribution was evaluated. The neutron (E_th<0.5 eV) distribution at the sample position was calculated. A mesh tally with pixels 2 mm×2 mm was mapped to calculate the distribution. The results are presented in Fig. 3(a).

Fig. 3Full-screen View

(Color online) (a) Mapping of thermal neutron flux. (b) Typical gamma ray spectra and fitting peaks. (c) Ratio of Experiment/MCNP for each region

Benchmark experiments were conducted to verify the simulation results. Owing to the complicated process and low efficiency of foil activation, an alternative method based on a Cd-Zn-Te detector was used to determine the thermal neutron distribution at the sample position. Because the induced gamma-rays are produced inside the crystal, the thermal neutron detection efficiency is relatively high [39]. A Cd-Zn-Te detector with dimensions of 1 cm×1 cm×0.5 cm (DT-01C1, Imdetek) was used to characterize the thermal neutron distribution. Considering symmetry, the measurements were conducted in only one-quarter of the region. As shown in Fig. 3(a), the measurement regions were divided by 3×3 parts (red box) based on the size of the detector.

The acquisition time of the Cd-Zn-Te detector was set to 600 s for each measurement. Typical gamma ray spectra and fitting peaks are shown in Fig. 3(b). The 558 keV peak was used for comparison with the simulated data. Net peak counts were obtained by subtracting the background signals. For the simulated data, the pixel values of each region were averaged and normalized to the highest experimental value. The results are presented in Fig. 3(c). We observed that the discrepancies were within 10%. A satisfactory agreement between the experimental and simulated data was observed, which indicates the workability of the MCNP simulation.

Model of HPGe detector

In this study, an N-type HPGe detector (ORTEC: TRANS-SPEC) with a relative efficiency of 55% and a resolution of 1.9 keV at 1333 keV was used to detect prompt gamma rays. The signals were collected using a multichannel analyzer and the MAESTRO software [40]. The structural parameters supplied by the detector manufacturer indicated that the diameter and length of the Ge crystal were 67 and 71 mm, respectively. There was an internal hole whose diameter and length are 15 and 60 mm, respectively. The thicknesses of dead layer of the internal hole and the external crystal surfaces are 700 and 0.3 μm, respectively. The Ge crystal was positioned inside an aluminum case with a diameter of 83 mm and length of 106.5 mm. As the parameters of the aluminum case were not provided by the manufacturer, an X-ray radiography system was used to observe the internal position of the crystal and the dimensions of the case. Figure 4(a) shows an image of the HPGe detector.

Fig. 4Full-screen View

(Color online) (a) X-ray radiography of HPGe detector. (b) MCNP model. (c), (d) Simulated and experimental gamma ray spectra of ¹³⁷Cs and ⁶⁰Co

An MCNP model of the HPGe detector head was constructed based on these parameters, as shown in Fig. 4(b). The simulated gamma ray spectrum was calculated using the pulse height tally (F8) and Gaussian energy broadening (GEB) cards. To verify the feasibility of the developed detector model, benchmark measurements were performed with reference to the ¹³⁷Cs and ⁶⁰Co point sources. The activities of the ¹³⁷Cs and ⁶⁰Co sources at the time of measurement were 9.19 × 10⁴ and 1.20 × 10⁴ Bq, respectively. The emission probabilities of 662 keV (¹³⁷Cs), 1173 keV and 1332 keV (⁶⁰Co) are 0.851, 0.9985 and 0.9998, respectively. The distance between the source and detector surfaces was 15 cm. The acquisition time was 900 s for the live-time model.

The experimental and simulation results are presented in Figs. 4(c) and (d). The simulation results agreed well with the experimental data. The full energy peak efficiencies (FEPEs) of the three peaks were compared. The discrepancies between the experimentally determined and simulated FEPEs were 3.76%, 3.71% and -2.18%, respectively. The main reasons are probably counting statistics and scattering by surrounding materials. The developed model was validated based on these uncertainties.

Sample measurement and image analyses

To ensure a satisfactory counting statistic, the measurement time was set to 3 h with a live time due to the low thermal neutron flux. The induced gamma ray spectra of were also recorded using the MAESTRO software. The analysis of the prompt gamma ray peaks was carried out using the GAMMAFIT software [41] because the interference and overlapping of the peaks caused a challenge for accurate peak fitting owing to the poor SNR. In addition, a traditional method with a single collimated detector was also performed; the measurement time was the same, and the square opening was also 1 cm×1 cm. After setting the measurement time, the neutron damage to the HPGe detector was evaluated by calculating the accumulated neutron fluence. An F4 tally was used to calculate the neutron flux in the HPGe crystal. The average value was approximately 4.71×10^-5 cm^-2 and was multiplied by the neutron yield and measurement time. The result was 7.36×10⁶ cm^-2, which is acceptable compared to published work [42, 43].

The image quality was evaluated using the structural similarity index measure (SSIM) and root mean square error (RMSE) [44]. SSIM is a widely used metric that reflects the visual similarity between a true image and a reconstructed image; an SSIM value close to unity indicates that the two images are similar. The SSIM is defined as follows: $\text{SSIM}=\frac{\left(2{\mu }_{\text{r}}{\mu }_{\text{t}}+C_1\right)\left(2{\sigma }_{\text{r,t}}+C_2\right)}{\left({\mu }_{\text{r}}^2+{\mu }_{\text{t}}^2+C_2\right)\left({\sigma }_{\text{r}}^2+{\sigma }_{\text{t}}^2+C_2\right)},$ (8) where μ_r, μ_t, σ_r and σ_t are the average and standard deviations of the reconstructed and true images, respectively; σ_r,t is the covariance between the reconstructed and true images; C₁ and C₂ are constants set to avoid zero denominators. In this study, C₁ and C₂ were set to 0.1.

RMSE reflects the difference in pixel values between a true image and a reconstructed image. A smaller value indicates that the reconstructed image has a smaller statistical deviation. The RMSE is defined as follows: $\text{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N}{\left(P_i^{\text{r}}-P_i^{\text{t}}\right)}^2}\mathrm{,}$ (9) where N is the pixel number, $P_i^{\text{r}}$ and $P_i^{\text{t}}$ are the ith pixel values of the reconstructed and true images, respectively.

Furthermore, the spatial resolution of the proposed system was calculated using the edge spread function (ESF), which describes the spreading or blurring of edges or sharp transitions in an image [45]. This quantifies how the intensity values change across an edge, providing information on the quality of image sharpness. The ESF is obtained by analyzing the pixel values along a line perpendicular to an edge or transition in an image. The line is usually taken across the transition from background to foreground, or vice versa. By analyzing the intensity values along this line, the ESF represents the changes in intensity from one side of the edge to the other.

Analysis of induced gamma ray spectra

The enlarged gamma ray spectra of the region of interest for the pattern collimator and the traditional method are shown in Fig. 5(a) and (b). The spectra essentially consist of induced gamma rays of elements from the sample and other materials, including the collimator, shielding, and detector materials. The major isotopes and their gamma ray energies are marked. The most prominent peaks were produced by iron and aluminum through the inelastic scattering reactions ⁵⁶Fe (n,n’γ) ⁵⁶Fe and ²⁷Al (n,n’γ) ²⁷Al. These gamma rays originated mainly from the collimator and detector housing materials, which correspond to the 847 keV and 1238 keV of the Fe peaks, as well as 844 keV and 1014 keV of the Al peaks. For the 962 keV of the Cu peak, they were mainly emitted from the detector cold finger. Moreover, Pb peaks could be observed because they were components of the shielding material. Another fraction of induced gamma rays was emitted from the detector crystal, the neutrons reacted with the germanium nucleus which led to the 608 keV, 868 keV and 1204 keV peaks through INC reaction, as well as the formation of the so-called germanium triangles 597 keV, 696 keV and 834 keV through INS reaction.

Fig. 5Full-screen View

(Color online) (a) Schematic diagram of two methods (b) Gamma ray spectra in the energy range 520-1250 keV acquired during 10800 s for this method (Orange) and the traditional method (Red). (c) Fitting of Cd peak and Cl peak. (d) Count rates and SNRs of Cl peak and Cd peak for each collimator

The overall count rates in the interesting region were 8160 cps for the proposed method and 1800 cps for the traditional method. This is because more square openings lead to more background gamma-rays and neutrons interacting with HPGe detector. Although the background signal was approximately four times higher than that of the traditional method, the peaks of Cd and Cl were still readily observed. The 558 keV and 748 keV Cd peaks were identified. However, the 651 keV Cd peak could not be observed because it was affected by the 597 keV germanium triangle. For the gamma rays emitted from Cl, the 788 keV, 1164 keV, 1171 keV and 1173 keV peaks were identified. All of these peaks could not be observed in the spectra acquired using the traditional method. As mentioned above, the intensities of the 558 keV Cd and 1164 keV Cl peaks were analyzed using the GAMMAFIT software. An example is shown in Fig. 5(c). The count rates and corresponding SNRs for the two peaks in the 30 measurements are shown in Fig. 5(d). Under the present measurement conditions, the average SNRs were approximately 8% for the 1164 keV Cl peak and 12% for the 558 keV Cd peak.

Correction of neutron field

For neutron-induced gamma ray imaging, the neutron field inside a sample should ideally be maintained. However, as shown in Fig. 3(a), the thermal neutron flux decreases at the edge owing to the isotropic emission of ²⁴¹Am-Be neutron source. Moreover, the neutron beam is attenuated, particularly for elements with high TNC macroscopic cross-sections, which leads to a neutron self-shielding effect. To solve this problem, neutron radiography/tomography can be used to experimentally determine the neutron field. However, neutron imaging using a low-yield source is not feasible. An alternative approach is to calculate the neutron flux inside each pixel using Monte Carlo simulation [46]. Correction of the solid angle and neutron self-shielding effect can be performed as follows: $P_i^{\text{C}}=P_i^{\text{r}}\frac{{\varphi }_{i\mathrm{,}\text{d}}}{{\varphi }_{i\mathrm{,}\text{r}}}\frac{{\varphi }_{\text{f,max}}}{{\varphi }_{i\mathrm{,}\text{f}}}$ (10) where $P_i^{\text{C}}$ and $P_i^{\text{r}}$ are the i_th pixel values of the corrected and reconstructed images, respectively, and ${\varphi }_{i\mathrm{,}\text{f}}$ and ϕ_f,max are the ith pixel and maximum thermal neutron flux shown in Fig. 3(a). ${\varphi }_{i\mathrm{,}\text{r}}$ is the thermal neutron flux inside the ith pixel, ${\varphi }_{i\mathrm{,}\text{d}}$ is the one for “diluted sample,” where the density of sample is reduced by a factor of 1000. Neutron beam attenuation is negligible and reflects neutron self-shielding [47]. These values are obtained using a mesh tally.

Image reconstruction

After obtaining the net peak areas of the 30 measurements, the Cd and Cl images were reconstructed using the MLEM algorithm. The SSIM and RMSE for different numbers of iterations were calculated, as shown in Fig. 6(a) and (b). For SSIM, the optimal number of iterations was 20 (Cl images) and 70 (Cd images). For RMSE, the optimal number of iterations was 40 (Cl image) and 70 (Cd images). This was primarily because the effect of noise on the reconstructed image tended to increase with the number of iterations when the SNR was small. Finally, the number of iterations is set to 70. The corresponding reconstructed images are presented in Figs. 6(c) and (d). The color bar indicates the normalized content, ranging from 0 (minimum value) to 1 (maximum value). The SSIM values of the Cd and Cl images are 0.67 and 0.65, respectively, whereas the RMSE of the Cd and Cl images are 0.36 and 0.37, respectively. For easier visualization, as shown in Fig. 6(e) shows the red-green-blue (RGB) map of the elemental distributions. The reconstructed images show the elemental distribution of the sample accurately. The results show that the distributions of Cd (red) and Cl (green) can be clearly distinguished. The regions occupied by Cd are observed in the four corners, whereas the remaining regions correspond to Cl.

Fig. 6Full-screen View

(Color online) (a) SSIM as a function of the iteration number. (b) RMSE as a function of the iteration number. (c), (d) Reconstructed element imaging of the sample. The reconstructions on the left and right are the Cd and Cl images, respectively. (e) RGB map of the sample, where red represents Cd and green represents Cl. (f) ESFs of Cd and Cl images

Figure 6(f) shows the ESFs and corresponding fitting functions of two element images. The ESF function is obtained by averaging the pixel values within the yellow box in Fig. 6(e), horizontally followed by fitting with a sigmoid curve (orange line). Subsequently, the first derivative of the sigmoid curve is calculated and fitted using a Gaussian function (blue line). The full width at half maximum (FWHM) of the Gaussian curve is defined as the spatial resolution. The results show that the spatial resolutions of Cd and Cl images are 0.66 cm and 0.62 cm, respectively. The results reported by Laszlo [48] and Kluge [49], obtained using a collimated neutron beam at the BNC and FRM II reactors, were approximately 0.2–0.3 cm. The result reported by Chen Mayer [23] using a Compton camera detector on the NIST reactor was approximately 0.3 cm. Compared with these studies, we achieved a similar level of resolution using a low-flux neutron source.

In this study, the element images were reconstructed using a relatively low flux, resulting in data with a low SNR and a long measurement time. In practical applications, the measurement time can be reduced by using a neutron source with a higher neutron yield. Scintillator detectors (NaI, BGO, LaBr₃, etc.) combined with a full gamma ray spectrum analytical algorithm can also be used in certain specific cases owing to their high detection efficiencies and counting statistics. In addition, the collimator masks can be optimized to improve the quality of the reconstructed element images when working at a low SNR or low sampling rates. A rotating collimator can be used because the modulation varies with the rotation angle. The incident gamma rays produce a fluctuating signal corresponding to the modulation at different angles [32].

When applying this technique to real objects with unknown information, it must be combined with other techniques for a comprehensive analysis. The acquisition of accurate geometrical information and chemical compositions is vital when constructing a Monte Carlo model. Optical scanning and X-ray tomography are effective techniques for obtaining accurate geometrical information [50]. An initial scan measurement can be conducted over a relatively large volume to provide a rough estimate of chemical composition. Based on the data, a preliminary Monte Carlo model was developed. Subsequently, its parameters were adjusted to accurately simulate the detection efficiency matrix and neutron field. This fine-tuning process involves iterations in which the compositions assigned to individual pixels are continuously refined. Through this iterative refinement, a state of self-consistency was achieved, thereby generating an elemental map with a high spatial resolution [47].