Optimization of transmission reconstruction equation

In traditional TGS transmission measurements, the transmission reconstruction equation is expressed as Eq. (1) [22]: ${\displaystyle \sum _{j=1}^n{\mu }_j}{T}_{ij}=p_i$ (1) where pi is the negative logarithm of the γ-rays transmittance at the i-th transmission measurement position, also known as projection. Further, Tij is the trajectory length of the j-th voxel that the γ-rays pass through at the i-th position, μj is the line attenuation coefficient of the j-th voxel, and n is the total number of voxels.

In actual TGS transmission measurements, γrays emitted by the transmission source exhibit an angular distribution, and the detection efficiencies of all points on the detector end face are not the same. The TGS transmission formula is expressed as follows [23]. $p_i=\frac{1}{{\displaystyle {\sum }_{k=1}^N{\epsilon }_k}}{\displaystyle \sum _{k=1}^N\left[{\epsilon }_k\exp \left({\displaystyle \sum _{j=1}^n\left(-x_{ikj}{\mu }_j\right)}\right)\right]}$ (2) where pi is the transmissivity at the i-th position, εk is the detection efficiency of the k-th point on the detector end face for a certain energy γ-rays, xijk is the length of the trajectory where the ray, passing through the k-th point of the detector end face, intersects the j-th voxel at the i-th measurement position, and N is the total number of γ-rays emitted by the transmission source within the solid angle of the detector.

Further, $q_k=\exp \left({\displaystyle \sum _{j=1}^n\left(-x_{ikj}{\mu }_j\right)}\right)$ (3) where qk is the attenuation rate at which the kth ray passes through the material without considering the detector efficiency.

Thus, Eq. (4) can be obtained as $p_i{\displaystyle \sum _{k=1}^N{\epsilon }_k}={\displaystyle \sum _{k=1}^N{\epsilon }_k}{q}_k$ (4) Next, we obtain Eq. (5) by considering the negative logarithms of both sides of Eq. (3): $-\ln q_k={\displaystyle \sum _{j=1}^n{x}_{ikj}}{\mu }_j$ (5) Then: $\left[\begin{array}{ccccc}{x}_{i11}& x_{i12}& x_{i13}& \dots & x_{i1n}\\ x_{i21}& x_{i22}& x_{i23}& \dots & x_{i2n}\\ ⋮& ⋮& ⋮& \dots & ⋮\\ ⋮& ⋮& ⋮& \dots & ⋮\\ x_{iN1}& x_{iN2}& x_{iN3}& \dots & x_{iNn}\end{array}\right]\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ ⋮\\ {\mu }_n\end{array}\right]=-\left[\begin{array}{c}\ln q_1\\ \ln q_2\\ ⋮\\ ⋮\\ \ln q_N\end{array}\right]$ (6) According to Taylor’s formula, Eq. (7) can be obtained as $\ln q_k=\left(q_k-1\right)-\frac{1}{2}{\left(q_k-1\right)}^2+\dots$ (7) Then: $\begin{array}{l}\left[\begin{array}{ccccc}{x}_{i11}& x_{i12}& x_{i13}& \dots & x_{i1n}\\ x_{i21}& x_{i22}& x_{i23}& \dots & x_{i2n}\\ ⋮& ⋮& ⋮& \dots & ⋮\\ ⋮& ⋮& ⋮& \dots & ⋮\\ x_{iN1}& x_{iN2}& x_{iN3}& \dots & x_{iNn}\end{array}\right]\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ ⋮\\ {\mu }_n\end{array}\right]\\ =-\left[\begin{array}{c}\ln q_1\\ \ln q_2\\ ⋮\\ ⋮\\ \ln q_N\end{array}\right]\approx -k\left[\begin{array}{c}{q}_1-1\\ q_2-1\\ ⋮\\ ⋮\\ q_N-1\end{array}\right]\end{array}$ (8) where $k=\frac{\ln \left(p_i/{\eta }_i\right)}{p_i/{\eta }_i-1}$ is the error correction value and ηi represents the detection efficiency of the detector for a certain energy γ-rays at the i-th measurement position.

However, Eq. (8) does not consider endpoint detection efficiency εk. If both sides are left multiplied by the detection efficiency [28], then we obtain $\begin{array}{l}\left[\begin{array}{cccc}{\epsilon }_1& {\epsilon }_2& \dots & {\epsilon }_N\end{array}\right]\left[\begin{array}{ccccc}{x}_{i11}& x_{i12}& x_{i13}& \dots & x_{i1n}\\ x_{i21}& x_{i22}& x_{i23}& \dots & x_{i2n}\\ ⋮& ⋮& ⋮& \dots & ⋮\\ ⋮& ⋮& ⋮& \dots & ⋮\\ x_{iN1}& x_{iN2}& x_{iN3}& \dots & x_{iNn}\end{array}\right]\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ ⋮\\ {\mu }_n\end{array}\right]\\ =-k\left[\begin{array}{cccc}{\epsilon }_1& {\epsilon }_2& \dots & {\epsilon }_N\end{array}\right]\left[\begin{array}{c}{q}_1-1\\ q_2-1\\ ⋮\\ ⋮\\ q_N-1\end{array}\right]\end{array}$ (9) Eq. (9) is transformed as follows: $\begin{array}{l}{\mu }_1{\displaystyle \sum _{k=1}^N{\epsilon }_k}{x}_{ik1}+{\mu }_2{\displaystyle \sum _{k=1}^N{\epsilon }_k}{x}_{ik2}+\dots +{\mu }_n{\displaystyle \sum _{k=1}^N{\epsilon }_k}{x}_{ikn}\\ =-k[\left({\epsilon }_1{q}_1-{\epsilon }_1\right)+\left({\epsilon }_2{q}_2-{\epsilon }_2\right)+\dots \\ +\left({\epsilon }_N{q}_N-{\epsilon }_N\right)]\end{array}$ (10) By substituting Eq. (4) into Eq. (10), we obtain $\begin{array}{l}{\mu }_1{\displaystyle \sum _{k=1}^N{\epsilon }_k}{x}_{ik1}+{\mu }_2{\displaystyle \sum _{k=1}^N{\epsilon }_k}{x}_{ik2}+\dots +{\mu }_n{\displaystyle \sum _{k=1}^N{\epsilon }_k}{x}_{ikn}\\ =k\left(1-p_i\right){\displaystyle \sum _{k=1}^N{\epsilon }_k}\end{array}$ (11) The equivalent trajectory length can be expressed as $T_{ij}=\frac{{\displaystyle {\sum }_{k=1}^N{\epsilon }_k}{x}_{ikj}}{{\displaystyle {\sum }_{k=1}^N{\epsilon }_k}}$ (12) The system of equations can be expressed as Eq. (13): ${\displaystyle \sum _{j=1}^n{\mu }_j}{T}_{ij}=\frac{\ln \left(p_i/{\eta }_i\right)}{p_i/{\eta }_i-1}\left(1-p_i\right)$ (13) The equation for the equivalent trajectory length, which considers the endpoint detection efficiency, is expressed as Eq. (12). By substituting the obtained equivalent trajectory lengths Tij, detection efficiency ηi, and transmissivity pi into Eq. (13), an equation system for the transmission image reconstruction can be obtained.

Model of virtual trajectory length based on Geant4

To solve the optimized transmission image reconstruction equation, the equivalent trajectory lengths, detection efficiency, and transmissivity must be obtained in advance. Transmissivity can be obtained through detector counting conversion. For detection efficiency, a passive efficiency calibration method based on the Monte Carlo method is utilized. This method requires the calculation of the detection efficiency of the detector relative to the radiation source for a certain energy gamma ray at each measurement position in advance based on the Monte Carlo method. However, studies on the solution of the equivalent trajectory length formula in Eq. (12) are scarce. The trajectory lengths obtained using Eq. (12) can be interpreted as the average length of all γ-rays detected by the detector passing through the jth voxel. If the trajectory lengths of all γ-rays passing through the j-th voxel can be calculated, then the fact whether each ray has been detected by the detector can be determined [24]. The equivalent trajectory length is obtained by averaging the trajectory lengths of the detected rays. The Monte Carlo program Geant4 can track the trajectory lengths of each interaction step between the primary or secondary particles and matter and obtain the geometric information of each step, which provides the possibility for simulations to obtain equivalent trajectory lengths [25, 26].

To determine the equivalent trajectory length and detection efficiency, this study proposed a novel approach by constructing a virtual waste barrel and integrating the solution processes into a unified model. This model, based on Geant4, simplified the solving process and introduced a pioneering method for determining equivalent trajectory lengths. The virtual trajectory length model was structured as follows:

Based on the TGS experimental setup, a TGS simulation model without a nuclear waste barrel was constructed. The space containing the virtual waste barrel was divided into multiple segments, and one segment was selected for the simulation to calculate the equivalent trajectory length [27]. The virtual waste barrel segment was further divided into N × N voxels, and the medium within the segment was filled with vacuum. The voxel segmentation of a single-layer nuclear waste barrel and emission status of the rays are shown in Fig. 1.

Fig. 1Full-screen View

(Color online) Schematic of the virtual trajectory length model based on Geant4. The model divided the segment into 14×14 voxels, with γ-rays energy set to Ek=1.17 MeV. This figure shows the situation wherein particles emitted from the radiation source passed through certain voxels. Based on Eq. (12), the model can obtain the number of particles detected by the detector on a certain voxels, as well as the cumulative sum of the trajectory lengths left by these particles on that voxel(${\displaystyle {\sum }_{k=1}^N{\epsilon }_k}{x}_{ikj}$)

The next step involved simulating the transportation of numerous particles passing through a virtual waste barrel segment at a specific location. Leveraging Geant4 functions such as “G4Step,” “Getname(),” “GetTouchableHandle(),” and “GetPostStepPoint(),” the name of voxels traversed by the particles and the length of trajectory within these voxels were obtained. Notably, because of the absence of any substances in the virtual waste barrel, the particle rays followed a straight-line transport path without undergoing any physical processes [29].

Considering the uneven response of the end- face detection efficiency of the detector, the simulation included a process wherein particles passing through the virtual waste barrel were detected by an HPGe detector. If the energy deposited in the detector differed from the particle’s initial energy, it was assumed that the particles were not detected, and the trajectory lengths passing through all voxels were considered as 0.

The cumulative trajectory lengths of all the particles with the same emission energy passing through the same voxel were then calculated. Dividing the total trajectory length by the number of particles with a nonzero trajectory length yielded an equivalent trajectory length corresponding to the voxel. Further, dividing the number of particles with a nonzero trajectory length at a specific energy by the total number of emitted particles at that energy yielded the detection efficiency of γ-rays at a particular location, denoted as ηi.

Image reconstruction algorithm

The current image reconstruction algorithms can be roughly divided into two categories: analytical methods and iterative reconstruction algorithms [30]. This study focused on accurately obtaining the transmission measurement trajectory length and did not delve deeply into the intricacies of the reconstruction algorithm. The optimization of the algorithm will be discussed in future research. Consequently, for image reconstruction, the Maximum Likelihood Expectation Maximization (MLEM) algorithm was employed [31], and its iterative formula is expressed as follows [32, 33]: ${\mu }_j^{\left(k+1\right)}=\frac{{\mu }_j^{\left(k\right)}}{{\displaystyle {\sum }_i{T}_{ij}}}{\displaystyle \sum _i{T}_{ij}}\frac{p_i}{{\displaystyle {\sum }_j{T}_{ij}}{\mu }_j^{\left(k\right)}}$ (14) where k is the number of iterations of the algorithm, ${\mu }_j^{\left(k\right)}$ is the estimated value of the line attenuation coefficients on all voxels during the k-th iteration, pi is the projection value obtained during the i-th transmission measurement of the substance to be measured, and Tij is the trajectory length value on the j-th voxel during the i-th measurement.

Image quality evaluation methods

To assess the quality of the reconstructed images of nuclear waste barrels, various evaluation methods are essential, including visual inspection and profiles for qualitative assessment, and the mean square error (MSE) and signal-to-noise ratio (SNR) for quantitative analysis [34].

MSE indicates the difference between the reconstructed and reference images, and a smaller MSE value indicates a better quality of the reconstructed image. In this study, SNR was also employed to evaluate the noise level of the reconstructed image; higher SNR values corresponded to lower image noise.

The MSE and SNR are expressed in Eq. (15), and Eq. (16): $MSE=\frac{1}{MN}{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{\left({\mu }_{\text{Rec}}\left(m\mathrm{,}n\right)-{\mu }_{\text{Ref}}\left(m\mathrm{,}n\right)\right)}^2}}$ (15) $SNR=10{{\displaystyle \log }}_{10}\left[\frac{{\displaystyle {\sum }_{m=1}^M{\displaystyle {\sum }_{n=1}^N{\mu }_{\text{Ref}}{\left(m\mathrm{,}n\right)}^2}}}{{{\displaystyle {\sum }_{m=1}^M{\displaystyle {\sum }_{n=1}^N\left({\mu }_{\text{Rec}}\left(m\mathrm{,}n\right)-{\mu }_{\text{Ref}}\left(m\mathrm{,}n\right)\right)}}}^2}\right]$ (16) where μ_Rec(m, n) is the reconstruction value of the mth row and nth column of the reconstructed image and μ_Ref(m, n) is the reference value of the mth row and nth column of the reference image.

Simulation model of nuclear waste barrel

To evaluate the virtual trajectory length model, Geant4 was used to simulate the TGS system. Two models were designed in this study. The first model, labeled Model 1, was a nuclear waste barrel without concrete filling. It comprised four distinct materials: fiber, concrete, water, and aluminum. The geometric representation of Model 1 is shown in Fig. 2(a). In addition, a second model, labeled Model 2, was developed to assess the practical applicability of the virtual trajectory length model. Model 2 was a nuclear waste barrel filled with concrete, as shown in Fig. 2(b). The fibers were positioned at three specific locations within Model 2, and the remaining space outside these regions was filled with concrete. Reference images of nuclear waste barrels are required to evaluate the quality of the reconstructed transmission images [35]. In this study, parallel-beam γrays were used to irradiate the nuclear waste barrel model vertically. The intensities of the rays passing through the model at different positions were obtained, which facilitated the acquisition of a reference image of the nuclear waste barrel [36]. In the reference image, different colors represent distinct densities of the medium [37], as shown in Fig. 2(c) and (d). Notably. in the projection calculations, the number of particles without attenuation is based on the number of particles passing through the barrel wall. The count after particle attenuation utilizes the particle count values of the particles passing through either Model 1 or Model 2. Consequently, the reconstructed image did not include the barrel wall, thus circumventing the potential influence of nuclear waste barrel walls on the results.

Fig. 2Full-screen View

(Color online) Geometric structure of Model 1 and Model 2, and their reference images

Comparison of different models for trajectory lengths

Two models were incorporated to validate the accuracy of the proposed trajectory length method. The first model was the “point-to-point” model established by Estep, and the second model was the average trajectory model proposed by Quanhu. These models provide additional trajectory length data for reference when a virtual model is used to extract trajectory lengths for image reconstruction. When reconstructing the images using trajectory lengths obtained from the “point-to-point” and average trajectory models, the projection data were not corrected for detection efficiency. However, when the virtual model was utilized for the trajectory length, the projection data were corrected for detection efficiency. The reconstruction images of Model 1 are labeled as Fig. 3(a), (c), and (e), and those of Model 2 are labeled as Fig. 3(b), (d), and (f).

Fig. 3Full-screen View

Reconstruction results of Models 1 and 2 under different trajectory length models. The results of Model 1 are shown in 1st column. The results of Model 2 are shown in 2nd column

Qualitative evaluation

In the reconstructed images (Fig. 3(a)) corresponding to the “point-to-point” model, the result appears unclear. They could not accurately capture the specific structure of the medium in Model 1. Conversely, the reconstruction image for the average-trajectory length model shown in Fig. 3(c) provides a distinct representation of the structure of the medium. Therefore, the average trajectory length model exhibited superiority over the “point-to-point” model. In the reconstruction images (Fig. 3(e)), corresponding to the virtual trajectory length model, a visual inspection confirmed that the reconstructed medium (both shape and position) within the barrel was closer to the reference image of Model 1. This indicated a substantial improvement in the reconstructed images.

In the reconstruction images of Model 2 (Fig. 3(b), (d), and (f)), it was evident apparent that, compared to the “point-to-point” and average trajectory length models, the virtual trajectory length model exhibited superior reconstruction quality. This model successfully reconstructed the filling material within the barrel and effectively restored pertinent information of the measured medium.

The reconstruction and reference images’ profiles of the two models are shown in Figs. 4(c), (e), and (d), (f), respectively. The positions of the profiles are presented in Fig. 4(a) and (b). A comparison of the images revealed that the reconstructed images’ profiles using the virtual trajectory length model exhibited a higher degree of consistency with the profile diagram of the reference image. The overall trends of their profiles were closer and smoother, implying that the proposed virtual trajectory length model significantly reduced the noise level of the reconstructed images. The reconstruction results of the virtual trajectory length model exhibited a qualitative improvement, whether for assumed Model 1 or Model 2, which was closer to the real waste barrel environment.

Fig. 4Full-screen View

(Color online) Comparison of reconstruction image’ profiles

Quantitative evaluation

For the various trajectory length models, the MSE and SNR of the reconstructed images for both Models 1 and 2 are shown in Fig. 5:

Fig. 5Full-screen View

(Color online) MSE and SNR of reconstructed images under different trajectory length models

In Fig. 5(a) and (b), the virtual trajectory length model demonstrated significant advantages over the “point-to-point” and average trajectory length models for Model 1. Compared to the “point-to-point” model, the virtual trajectory length model reduced the MSE by 57.8% and increased the SNR by 192.3%. Further, the average trajectory length model, which addressed the overlooked solid angle problem inherent in the “point-to-point” model, yielded trajectory lengths that were closer to the actual values. However, compared with the average model, the virtual trajectory length model proposed in this study still achieved a 37.6% reduction in MSE and a 56.1% improvement in SNR. This discrepancy can be attributed to the following factors. (1) The number of detector end-face grids. The process of numerous particles reaching the detector end face and being detected was simulated using a virtual trajectory length model. Notably, the number of end-face grids was considerably greater than that of the average model, reaching millions or even billions of levels (as determined by the number of emitted particles, more particles will lead to more grids). (2) Non-uniform detection efficiency. The proposed model considered the non-uniform detection efficiency of the actual TGS detector end face. This resulted in a smaller deviation between the obtained and actual trajectory lengths.

In the objective assessment of the reconstruction results for Model 2, the virtual trajectory length model consistently demonstrated a superior reconstruction quality. Compared to the “point-to-point” and average trajectory models, the MSE decreased by 85.4% and 72.5%, respectively. In addition, the SNR increased by 375.0% and 112.7%, respectively. These results confirmed that the trajectory length model proposed in this study improved the quality of the reconstructed images and effectively reduced noise.

Based on the comprehensive qualitative and quantitative analyses conducted above, it can be concluded that the proposed virtual trajectory length model, which considers the nonuniform detection efficiency of the detector end face and the angular distribution of γ-rays at the detector, could accurately reconstruct the distribution of the medium within the nuclear waste barrel. The reconstruction quality achieved by this model was significantly superior to that of both the “point-to-point” and average trajectory models.