Quantitative imaging and automated fuel pin identification for passive gamma emission tomography

Compliance of member States to the Treaty on the Non-Proliferation of Nuclear Weapons is monitored through nuclear safeguards. The Passive Gamma Emission Tomography (PGET) system is a novel instrument developed within the framework of the International Atomic Energy Agency (IAEA) project JNT 1510, which included the European Commission, Finland, Hungary and Sweden. The PGET is used for the verification of spent nuclear fuel stored in water pools. Advanced image reconstruction techniques are crucial for obtaining high-quality cross-sectional images of the spent-fuel bundle to allow inspectors of the IAEA to monitor nuclear material and promptly identify its diversion. In this work, we have developed a software suite to accurately reconstruct the spent-fuel cross sectional image, automatically identify present fuel rods, and estimate their activity. Unique image reconstruction challenges are posed by the measurement of spent fuel, due to its high activity and the self-attenuation. While the former is mitigated by detector physical collimation, we implemented a linear forward model to model the detector responses to the fuel rods inside the PGET, to account for the latter. The image reconstruction is performed by solving a regularized linear inverse problem using the fast-iterative shrinkage-thresholding algorithm. We have also implemented the traditional filtered back projection (FBP) method based on the inverse Radon transform for comparison and applied both methods to reconstruct images of simulated mockup fuel assemblies. Higher image resolution and fewer reconstruction artifacts were obtained with the inverse-problem approach, with the mean-square-error reduced by 50%, and the structural-similarity improved by 200%. We then used a convolutional neural network (CNN) to automatically identify the bundle type and extract the pin locations from the images; the estimated activity levels finally being compared with the ground truth. The proposed computational methods accurately estimated the activity levels of the present pins, with an associated uncertainty of approximately 5%.

to increase the sensitivity to a variety of fuel diversion scenarios 5 . The software to analyze PGET images can be improved 7 to reduce user intervention and overall perform faster and more reliable monitoring of spent fuel and therefore implement the ultimate safeguards objective, i.e., the prompt detection of nuclear material diverted from peaceful uses.
The PGET of spent fuel poses some unique challenges, compared to industrial or medical tomographic imaging, which are due to the high activity of the sources (a single-pin activity is of the order of 10 13 Bq) and the high self-attenuation of the fuel pins. While the former can be mitigated using collimated detectors, the latter needs to be addressed by using specific imaging algorithms. Previous work has shown that the reconstruction quality could be significantly improved by applying the attenuation correction 8 . In this work, we further corrected for the gamma ray down-scattering in the energy window of interest. We have implemented and integrated computational methods based on proximal gradient, physics-informed Monte Carlo sampling, and machine learning to (1) improve the PGET image quality, (2) automatically identify missing pins, and (3) quantify the pin activities. As we will show, the automated identification of missing pins relies heavily on the image quality, and therefore the three tasks are inherently intertwined.

Results
Simulated sinograms. We simulated the measurement of six fuel pin configurations in a mock-up waterwater energetic reactor (VVER) fuel assembly using MCNP 6.2 9 , as shown in Fig. 1. In Cases 2 and 3, we mimicked the scenario where fuel pins were missing. In Cases 4 and 5, we used depleted uranium to replace cobalt because it resulted in the highest attenuation among different potential replacement materials 10 , which would allow us to examine the algorithm's capability of distinguishing replaced pins. In Case 6, we wanted to explore whether the algorithm is able to render an accurate image of a space-dependent activity distribution. Details of the simulated PGET model are included in "Methods-Simulation Methods" section.
The detection system we simulated is the so-called PGET prototype device as described by 4 . The detection unit consisted of two collimated CdTe detector arrays on opposite sides of the fuel assembly, each encompassing 91 detectors. The detector arrays rotated and scanned the fuel bundle in steps of 1 • , which generated a 182 × 360 sinogram, as shown in Fig. 2. We simulated 5 × 10 9 NPS (number of particle histories) in each case and detected the photons in the 700-1500 keV energy window. Approximately 15,000 maximum counts in one pixel were achieved, comparable to the counts in an actual PGET measurement.
Inverse approach. We have developed a systematic approach that allows us to reconstruct high-quality images from individual sinograms, identify the pin locations, and quantify the pin activity levels. The image of fuel pins was reconstructed by solving a linear inverse problem. The region of interest was divided into 182 × 182 Figure 1. Actual pin distributions in the simulated fuel assemblies (ground truth). Pin activity is proportional to the brightness. In Case 1, the mock-up fuel assembly hosted 331 60 Co pins of 8.879 g/cm 3 density, which emitted 1.17 and 1.33 MeV gamma rays. In Cases 2 and 3, 10% 60 Co pins pins were removed at the center or uniformly. In Cases 4 and 5, 10% 60 Co pins pins were replaced by depleted uranium pins of the same size but higher density (10.4 g/cm 3 ). In Case 6, we decreased the activities of the middle pins by 80% and bottom pins by 50%.  www.nature.com/scientificreports/ pixels. We calculated the detector array response to a source of unit activity inside each pixel, forming the system response matrix. For an unknown fuel assembly, the measured sinogram is a linear combination of the calculated responses, with the coefficients being the source strength of each pixel. In this way, the image reconstruction problem was converted into a linear inverse problem, with both the simulated sinogram and calculated response matrix as inputs. First, we performed a simple image reconstruction to obtain an initial estimate of pin locations and pin radius. We calculated the system response matrix using a deterministic ray-tracing method, assuming that the scattering of gamma rays can be neglected 11 . The reconstructed image of Case 1 is shown in Fig. 3a. Inner pixels are brighter than the outer ones on the image, indicating an overestimation of the activity of the pins. This is because the contribution of scattered photons to the system response cannot be neglected for inner pixels, due to the heavy attenuation of unscattered photons. Nevertheless, using this image, we can determine the center of all possible pins as the centroids of the bright regions on the reconstructed image. In this step, we would allow some slackness in pin identification, since a more accurate reconstruction will be performed later. The identified pins are shown in Fig. 3b. Horizontal and vertical line profiles passing through the center of each pin were extracted from the image, and we fitted a Gaussian to the average of all profiles, shown in Fig. 4. The FWHM (full width at half maximum) of the fitted Gaussian was 0.769 cm, which was taken to be the pin diameter.
With the pin locations and pin radius known, we could create a material map of the fuel assembly. We then implemented an accelerated Monte-Carlo algorithm to calculate the system response matrix, accounting for both the absorption and scattering interaction of photons in the fuel assembly. The images reconstructed using the new response matrix are shown in Fig. 5. The images were then fed to a convolutional neural network to perform pin identification. The results are shown in Fig. 6. For Cases 1 to 5, we achieved 100% classification accuracy, due to the high quality of reconstructed images. We created a histogram of pin activities, as shown in Fig. 7. In Cases 1-5, the standard deviation of activity estimation ranged from 3 to 6%, and the mean of the pin activity distribution deviated from the ground truth by less than 4%. In Case 6, we successfully identified all pins, except those with the lowest activity level (Group 1) shown in Fig. 6f and Fig. 7f. The activity of pins in Group 2 and Group 3 was accurately reproduced with a relative error below 1%. The inverse approach is detailed in the "Methods-image reconstruction methods" section.

Discussion
The image reconstruction problem can be formulated into different linear inverse problems based on the selected model of observation noise. As detailed in "Inverse Problem Approach", we formulated two inverse problems with the Gaussian noise model and Poisson model, which can be solved using FISTA (fast-iterative shrinkagethresholding algorithm) 12 and PIDAL (Poisson image deconvolution by augmented Lagrangian) algorithm 13 ,  www.nature.com/scientificreports/ respectively. We applied both FISTA and PIDAL to perform reconstruction in Case 1 and compared them in Fig. 8. We assessed the image quality quantitatively using the mean-square-error (MSE) and structural similarity (SSIM) 14 . MSE is the mean-squared-difference between the reconstructed image and the ground truth, and SSIM measures the similarity between the reconstructed image and the ground truth. Lower MSE and higher SSIM mean better reconstruction. As shown in Table 1, the image reconstructed with FISTA resulted in a lower MSE and higher SSIM, compared to PIDAL, and PIDAL took approximately 20% longer to run than FISTA. Therefore, we used FISTA as the main image reconstruction algorithm. For comparison, we have also implemented the traditional FBP method to reconstruct the image, which is detailed in "Methods-image reconstruction methods" section. We then input these images to the neural network for pin identification, and estimated the pin activity levels. We compared the performance of the inverse (FISTA) and FBP approaches in terms of image quality, accuracy of pin identification and activity quantification. Fig. 9 Figure 5. Image reconstructed using the inverse approach. Good contrast between the pin-present region and pin-absent region was achieved. No severe artifacts were observed on the reconstructed images.   www.nature.com/scientificreports/ shows the images reconstructed using FBP. Compared to Fig. 9, images in Fig. 5 demonstrated higher contrast and fewer reconstruction artifacts. The blurring at the center was significantly removed using the inverse approach. We calculated the MSE and SSIM for both sets of images, shown in Table 2. We achieved significantly lower MSE and higher SSIM by using the inverse approach, which resulted in lower mis-classification rates of pin identification overall. The pin identification results based on FBP reconstruction are shown in Fig. 10. In Fig. 10a,c, and e, the blurring led to inaccurate pin localization around the center of FBP-reconstructed images, and the hexagonal structure was destroyed. In contrast, the hexagonal structure was accurately reproduced using the inverse problem approach, as shown in Fig. 6a-e. Hence, given the same sinogram data, the inverse approach is expected to better distinguish pins missing from the fuel assembly and result in lower false alarm rates, compared to FBP. MSE and SSIM are calculated based on unnormalized images. Therefore, in Case 6, we achieved the best MSE and SSIM because Case 6 exhibited the smallest average pixel intensity. Despite the good    www.nature.com/scientificreports/ image reconstruction, both FBP and inverse approaches led to relatively high mis-classification rates in Case 6 due to the fact that the low activity pins and high gradient in the pin activity are absent features in the training set. However, we do not expect such a large variation of pin activities in an actual fuel assembly 15 and, should that be the case, a new training set can be used to improve the pin identification results. We estimated the activity of each pin by summing the pixels on the reconstructed images and created a histogram of pin activities for each case, shown in Figs. 7 and 11. It should be noted that for the FBP reconstruction, there is no definite relationship between the pixel sum and the absolute activity. Appropriate normalization constant needs to be applied to convert the pixel sum into activity, which is not always possible in an actual inspection. In this case, we calculated the ratio between the mean of the pixel sum distribution and the true activity for Cases 2-5. We then used the average of these ratios as the normalization constant and applied it to all cases to convert the pixel sum to absolute activity. In contrast, normalization is not needed for the inverse reconstruction because the response matrix has already been normalized per unit source particle. Compared to Fig. 11, the pin activity distributions in Fig. 7 were narrower. We calculated the relative error compared to the ground truth and the standard deviation of pin activity, as shown in Table 3. We obtained smaller standard deviations by using the inverse approach in all cases. The relative error was negligible compared to the ground truth, except for the low activity pins in Group 1 of Case 6.
Compared to UO 2 fuel rod, the cobalt rod simulated in this work resulted in less attenuation because of its lower density and atomic number. Nevertheless, this work describes a general approach for correcting the gamma-ray attenuation and down-scattering in the energy window of interest and it can be easily adapted to UO 2 fuel assemblies by updating the atomic composition and gamma-ray source term in the simulation. In this work, we assumed that the detector performance is the same in the simulation of the sinogram and response matrix, which is not true for the real PGET device. The correction for varying detector performances across the detectors, either at the hardware or software level, is crucial to obtain an accurate estimation of pin activity. When identifying drifting detectors in post-processing, e.g., whose response is anomalously different compared to neighbour detectors, a simple approach consists in neglecting their response. In a preliminary analysis on simulated data, we replaced the response of an increasing number of detectors in the sinogram with null arrays. We found that if the number of "neglected" detectors is below ten in the PGET setup, the reconstructed image is negligibly affected by not including their response. Future work will be needed to study the robustness of the software suite as a function of a systematic bias in the detector performance.

Conclusion
In this work, we have implemented a full set of software, which is able to reconstruct cross-sectional images of mockup fuel assemblies acquired by a simulated PGET system, identify missing fuel pins, and estimate fuel pin activities based on the reconstructed image. We have developed a linear forward model that accounts for the scattering of gamma rays in the assembly to accurately characterize the response matrix of the PGET system. Figure 11. Activity distribution based on FBP. The pin activity was estimated by multiplying a normalization constant to the sum of all pixels inside the pin.The ground truth is shown as the red line. Table 3. Comparison of activity estimations. The mean NPS per pin is the mean of the activity distribution; relative error is the relative difference between the mean activity and ground truth; standard deviation refers to the standard deviation of the pin activity distribution. For FBP in Group 1 of Case 6, no data is shown because no present pin is identified. www.nature.com/scientificreports/ The image reconstruction was formulated into a linear inverse problem by modeling the observation noise as Gaussian, which was solved using FISTA. The reconstructed image was fed to a convolutional neural network to automatically identify the present pins and determine their centroids. Compared to the FBP approach, the inverse problem approach resulted in over 50% lower MSE and 200% higher SSIM, and consequently lower misclassification rates in pin identification in all cases. Based on the pin identification results, we estimated the pin activity by summing up the pixel values around the centroid inside the pin radius on the image. Compared to the FBP, the inverse approach resulted in smaller standard deviations of pin activity in all cases, with negligible bias with respect to the ground truth. The proposed inverse approach to reconstruct a fuel pin cross section, identify fuel pins and calculate their activity took approximately 8 minutes to run on an Intel Core i9-7920X CPU without parallelization. We are currently improving the algorithm to allow automatic classification of different activity groups in a real fuel assembly.

Methods
Monte Carlo simulation of the PGET based on MCNP. Figure 12 shows the cross-sectional view of the MCNP model of Case 1. In Case 4 and 5, depleted uranium pins were used, where Co was replaced by 0.20% enrichment UO 2 . We used the F4 tally as the detector response model in the MCNP simulation and the response matrix calculation. The F4 tally estimates the average photon flux in the detector cell by summing the track lengths of all particles in the 700-1500 keV energy window. The absolute activity measurement can be obtained using the proposed method as long as the simulated quantity in the response matrix is consistent with the simulated response to an unknown inspected fuel bundle. When applying the proposed technique to experimental data, one would want to validate the simulated detector response with the measured one, to properly incorporate specific detector's properties, such as its energy resolution and time response, in the simulated model.

Image reconstruction methods.
In this section, we describe in detail the two image reconstruction methods implemented in this work: the FBP and linear inverse problem method. We applied both methods to reconstruct images from the sinograms and compared their performance.
Filtered back-projection. Figure 13 shows the tomography measurement of the fuel assembly at angle θ . The FBP method relies on two approximations: first, we neglect the attenuation and scattering of gamma rays in the system; second, for pins at different locations, we assume that the geometric efficiencies are the same. Under these approximations, the sinogram φ(ρ, θ) can be simplified as the integral of the source distribution s(x, y) along the red line passing through the detector center, i.e., Equation (1) is the standard forward model in X-ray tomography imaging. The classical inverse operation from the sinogram to the source is the filtered back-projection, which is also a standard imaging algorithm in X-ray imaging 16 and will not be detailed here.
(1) φ(ρ, θ) = R 2 s(x, y)δ(x cos θ + y sin θ − ρ)dxdy. where A is the system response matrix of size M × N , and n models random observation noise, assumed to be isotropic Gaussian distributed. Physically, the i-th column of the response matrix A is the vectorized sinogram corresponding to a pin distribution with unit activity in pixel i but zero elsewhere. Accurate determination of the system response matrix is crucial to obtain a high-quality image and avoid systematic bias in pin identification and activity estimation. The analytical derivation of the response matrix is discussed in the next section. Given the simulated sinogram and system response matrix A , we can reconstruct the image by solving the following equation: where the first term is the data-fidelity term assuming Gaussian noise, and the second term is a regularization term acknowledging that the fuel pin is sparsely distributed. The regularization parameter is chosen based on the noise level. To solve Eq. (3), we have implemented the fast iterative shrinkage-thresholding algorithm (FISTA) 12 .
As an alternative, we can model the observation noise by Poisson noise and the data-fidelity term is changed accordingly 13,17 as follows To solve Eq. (4), we used the PIDAL algorithmic structure described in 13 .
Computation of the response matrix. The calculation of the system response matrix A involves the calculation of detector response to each pixel, with each response being one column in the matrix. For inner pixels, the contribution from scattered photons is non-negligible due to the high attenuation. To account for this, we have implemented an algorithm that combines stochastic photon transport with deterministic collimator-detector modeling to accurately calculate the system response matrix.
For real fuel assemblies, no prior information of pin distribution and pin radius will be given. It is therefore necessary to first estimate these quantities before running the Monte Carlo simulation. As discussed in the "Results-Inverse approach" section, this is done by performing a rough image reconstruction using our previously developed model, where uniform attenuation map is used and no scattering effect is considered 11 . Based on the reconstructed image, one could determine the pin locations and pin radius, and create a material map, where the material is assumed to be cobalt inside the pin, and water outside.
Next, we simulate the travel of photons inside the fuel assembly using Monte Carlo. Two sets of grids are generated. The fine grid has a pixel size of 0.5 × 0.5 mm 2 , while the coarse grid has a pixel size of 2 × 2 mm 2 . We first discretize the region of interest on the fine grid and calculate the response to the pixels whose center points are inside the pin. Then we calculate the response to a pixel on the coarse grid by summing the responses to the corresponding 4 × 4 pixels on the fine grid. Based on the pin identification result in the previous step, the potential pin-present region are pixelated into 5329 coarse pixels. We iterate over all these pixels and form the response matrix. In this way, we are able to reduce the pixelization error introduced during the discretization of the pins, without increasing the final response matrix size.
(2) φ = As + n, www.nature.com/scientificreports/ The response to each fine pixel whose center point is inside the pin is calculated in the following way. A photon (r 0 , � 0 , E 0 , W 0 ) is created, with the initial position r 0 uniformly sampled inside the pixel, the initial moving direction 0 uniformly sampled in 4π space, the initial energy E 0 sampled from the source energy spectrum, and the initial weight W 0 assigned to 1. We track the photon inside the PGET using the delta-tracking algorithm 18 and determine the position where the photon interacts with medium. We force the interaction to be Compton scattering by multiplying its weight with the probability that the interaction is Compton scattering, where µ sc (r i−1 , E i−1 ) and µ(r i−1 , E i−1 ) are the Compton scattering attenuation coefficient and total attenuation coefficient at the i-th interaction site, respectively. A scattering angle is sampled 19 , and the photon energy E i and moving direction i are updated accordingly. This process is repeated until the photon escapes the fuel assembly or the photon energy is below the detection threshold. A copy of the photon is saved at its creation site and each interaction site to a sub-projection map.
Once the sub-projection map is obtained, we apply the convolutional forced detection (CFD) technique to calculate the detector response. CFD is a widely-used variance reduction technique for X-ray downscatter simulation in single photon emission computed tomography (SPECT), which has been shown to be 50-100 times faster than conventional forced detection technique 20 and several thousand times faster than full-3D Monte Carlo simulation 21 . Figure 14 illustrates its principle. In CFD, the photon is forced to travel perpendicularly to the detector plane upon its creation or interaction with medium to effectively simulate all the possible spatial distributions of particle interactions and consequent detection events. For an unscattered photon, the contribution to the F4 tally in the detector cell is given by and for a scattered photon, the contribution is In Eqs. (6) and (7), W is the photon weight, is the solid angle subtended by the detector in steradian, σ (θ) is the differential cross-section for Compton scattering, µ(r, E) is the total attenuation coefficient along the CFD path, and the last term stands for the contribution to F4 tally, i.e., average track length in detector cell normalized by detector volume.
Due to the non-ideal collimation, a photon will contribute counts not only to the detector in the same column, but also to the neighbouring ones. This phenomenon is characterized by the depth-dependent point spread function (PSF) of the collimator, i.e., the photon count distribution on the detection plane when a point source is imaged. PSF depends on the vertical distance z between the photon and the detector array 22 : where b is the collimator aperture, l is the collimator length, µ is the total attenuation coefficient in the collimator, and FWHM is the width of the collimator PSF.
To calculate the response to a pin at an angle θ , the photons saved in the sub-projection map are first rotated by an angle −θ . We then sum the weights of the photons in the same pixel on the sub-projection map and calculate their contribution to the detector response based on Eqs. (6) or (7), which creates a projection map. We then  Figure 14. Illustration of the CFD technique. A photon shown as the gray disk is created on the left and travels in the medium along the solid arrow. The disk radius represents the photon energy. Upon it creation and each interaction with the medium, the photon is forced to be scattered along the dash arrow and to be detected by the detector with a certain probability. www.nature.com/scientificreports/ convolve each row in the projection map with the PSF function and sum all rows up to get the response at this angle. Next, we increment the angle θ and calculate the response at the next angle. In this way, photon tracking in the fuel assembly is done only once, which saves computation time.
We simulated the travel of 6400 photons to create the sub-projection map for one pixel. It took approximately 4 mins to calculate the response matrix containing 5329 pixels using CFD on an Intel Core i9-7920X CPU.
Pin localization. Images reconstructed using both approaches are input to a convolution neural network, referred to as U-net 23 , to perform pin identification. The experimental data provided by IAEA within the framework of a technology open challenge 7 were used for the training of the U-net. The training data consisted of the FBP images reconstructed from the measured sinograms of mock-up assemblies of 12 different geometries with a variable number of 60 Co pins, and examples of the reconstructed images of these assemblies can be found in 7 . The U-net was trained to minimize a binary crossentropy loss between the predicted pin masks and ground truth. For each image, we provided four randomly rotated ones as input to the network, for data augmentation. Data augmentation allowed us to generate new training instances from existing ones, artificially boosting the size of the training set. The network was trained for 10,000 epochs, i.e., complete iteration over the whole dataset, using a learning rate of 10 −5 and an Adam optimizer 24 . We found that 10 −5 is the learning rate for which the slope of the crossentropy loss is maximized, allowing to reach the minimum. The training was early stopped when the loss reached its minimum value. The trained U-net can extract visible structures from the input image and output a mask based on the extracted features. Bright regions on the mask with areas above a certain threshold are identified as present pins, and their centroids are calculated as the coordinates of the center of each pin.