Free-moving Quantitative Gamma-ray Imaging

The ability to map and estimate the activity of radiological source distributions in unknown three-dimensional environments has applications in the prevention and response to radiological accidents or threats as well as the enforcement and verification of international nuclear non-proliferation agreements. Such a capability requires well-characterized detector response functions, accurate time-dependent detector position and orientation data, a digitized representation of the surrounding 3D environment, and appropriate image reconstruction and uncertainty quantification methods. We have previously demonstrated 3D mapping of gamma-ray emitters with free-moving detector systems on a relative intensity scale using a technique called Scene Data Fusion (SDF). Here we characterize the detector response of a multi-element gamma-ray imaging system using experimentally benchmarked Monte Carlo simulations and perform 3D mapping on an absolute intensity scale. We present experimental reconstruction results from hand-carried and airborne measurements with point-like and distributed sources in known configurations, demonstrating quantitative SDF in complex 3D environments.

The reconstruction of the absolute intensity distribution of radioactive materials in complex 3D scenes has numerous applications including geologic survey, contamination assessment, guidance for decontamination and remediation verification, nuclear facility decommissioning, radiation safety/protection, and nuclear security. For example, 3D distributed source mapping is critical in response to radiological contamination such as in the Fukushima Prefecture in Japan, Chernobyl in Ukraine, or more recently the University of Washington Harborview Training and Research Facility where a sealed 2.9 kCi 137 Cs source was inadvertently breached 1 . In these cases, the assessment of the strength and distribution of the contamination was and remains to be necessary to minimize operator risk as well as to enable more effective means in the communication and visualization of nuclear radiation with the public to inform evacuation, decontamination, resettlement, and reopening efforts 2 . Additionally, such a capability can aid in consequence management efforts in order to avoid exposure to or plan dose-minimizing paths through radioactive fallout 3 . Furthermore, this capability facilitates search applications of compact radiological sources such as lost sources 4 , broken arrows (i.e., accidents involving nuclear weapons and components) 5 , as well as specific missions of the U.S. Special Operations Command such as sensitive site exploitation and nuclear disablement 3,6 .
The work presented here is part of an effort undertaken at Lawrence Berkeley National Laboratory (LBNL) and the University of California, Berkeley to develop a robust framework for performing 3D radiological imaging using platform-agnostic and free-moving detection and imaging systems. The emphasis of this paper is placed on quantitative gamma-ray image reconstruction, in that 3D gamma-ray maps are reconstructed accurately with relevant units such as activity of an isotope of interest (in µ Ci or MBq). In principle, the activity maps can also be used to compute 3D dose-rate maps (in µSv · hr −1 ) relevant for contamination assessment, emergency response, and radiation protection. Quantitative assessments are performed across scenes and environments (i.e., not just at the location of the instrument) and systems can be deployed on unmanned and remotely operating platforms to minimize risk to human operators if necessary. This paper begins by presenting the necessary background and overview of the quantitative free-moving 3D imaging approach, including data fusion methods and detector response characterization. Experimental demonstrations for both point-like and distributed gamma-ray sources are described and the results are discussed -highlighting the successes, current limitations, and potential areas for future work. Details regarding the gamma-ray imaging system, detector response characterization, image reconstruction algorithms, and uncertainty quantification are provided in the "Methods" section. www.nature.com/scientificreports/ Quantitative reconstruction Prior work with the SDF approach focused primarily on the development of the free-moving 3D imaging methods for various contextual and radiation sensor combinations and imaging modalities. The detector response functions used in reconstructions were not simulated or measured with respect to an absolute intensity scale and thus the gamma-ray activity estimates associated with the reconstructed images were qualitative in nature. The units associated with a gamma-ray image are critical to provide actionable information to operational personnel in applications such as contamination assessment, remediation verification and nuclear safeguards. Additionally, absolute units are necessary in biomedical imaging to provide quantitative assessment of disease and/or effectiveness of radiation-based treatments and therapies. In this work, we characterize the detector response functions with experimentally benchmarked Geant4 [31][32][33] simulations in order to produce quantitative images (see "Methods").
Experimental demonstrations Two source scenarios are considered and presented here. As an example of a port security application, the first scenario demonstrates the localization and quantification of a single pointsource in a wide-area outdoor environment involving stacks of cargo containers, using a small Unmanned Aerial System (sUAS)-based operation to move the radiation detection payload. The second scenario demonstrates the reconstruction of multiple distributed sources in a cluttered indoor laboratory environment, modeling a radioactive spill that can be surveyed using a detector system in a hand-held operation.
The SDF detector platform used in this work is MiniPRISM 34 -a 6 × 6 × 4 array partially populated with 1-cm 3 coplanar-grid CdZnTe (CZT) detectors 35,36 combined with a contextual sensor suite and on-board computer to perform lidar-based SLAM with Google Cartographer 37 . At the time of this measurement, a total of 58 detectors were loaded in the array in an optimized configuration to balance detection sensitivity and active coded mask imaging performance 38 . The MiniPRISM detector system and a rendering of the simulation model are shown in Fig. 1.

Results
Point source localization. A 1.84 mCi 137 Cs point-source was placed ≈2.5 m above the ground on the exterior of a steel cargo container in a mock port setting with >20 containers, some stacked two containers high (see Fig. 2a). The uncertainty of the source activity was undocumented, though we assume a 5-10% uncertainty (where uncertainty here is at 95% confidence or 2 σ ) based on previous purchases from the manufacturer. The MiniPRISM system was mounted to a sUAS and flown around the containers, mapping an area of >4,000 m 2 in under 7 min. The distance of closest approach to the source was ≈ 3 m. The average and maximum speeds of the sUAS during the measurement were 1.4 and 4.6 m · s −1 , respectively.
In contrast to online operation where multiple reconstructions can be performed as data are collected in realtime, the results presented here were generated as a single reconstruction that includes all data over the course of the measurement. The reconstruction energy region-of-interest (ROI) was set to 662 ± 60 keV to encompass the entire 137 Cs photo-peak for both singles (active coded mask) and doubles (Compton) events. While the expected energy resolution of coplanar-grid CZT is 1-3% 39 , the large energy ROI was necessary due to nonlinearities and gain drifts across individual detectors during the measurement. Here the doubles events refer to two-interaction sequences within the detector system in a coincident time window of 36 µ s. A total of 2864 singles and 370 doubles were recorded in the ROI over the entire measurement. Other reconstruction parameters include the following: 25 cm voxels, occupancy threshold of 5 points from the point cloud, 97,310 occupied voxels (1.5% of total), 1,916 poses (separated in time at 5 Hz) and 20 PSL iterations. At the time of the measurement, 7 of the 58 detectors were suffering from electronics problems and the counts were excluded from the reconstruction (though the occlusion provided by the detector material was still accounted for in the detector response).
The number of iterations used in maximum likelihood problems typically represents a balance between contrast recovery and image noise amplification 40 . Statistical stopping criteria have been explored, but generally do not translate to image quality metrics in terms of noise or artifacts 41 . Currently the number of iterations is set based on past experience and reconstruction speed requirements.   Fig. 2b where the localization estimates are displayed as Gaussian z-scores that represent the result of a likelihood ratio test between each voxel and the bestfit voxel (marked with a red star). The voxelized image estimate is overlaid on the measured point cloud, colorized  www.nature.com/scientificreports/ by the lidar return intensities, and the system trajectory, colorized by the gross counts in the detector. A topdown projection of the entire measurement space is shown, as well as three projections zoomed in near the true source location. Larger z-scores represent a larger deviation from the best fit and represent the spatial uncertainty of the best-fit estimate. In this visualization, z-scores above 5 σ are clipped from the image. The true source location (marked with a green dot) is 30 cm from the center of the best-fit voxel and is contained within the 2 σ spatial uncertainty. The estimated activity is 1.82 +0.13 −0.13 (stat.) +0.17 −0.17 (sys.) mCi of 137 Cs (2σ uncertainty bounds), and contains the true source activity of 1.84 mCi. The statistical uncertainty is calculated in the PSL algorithm (see "Methods") and the systematic uncertainty is calculated from a 9.4% systematic uncertainty (2σ ) ascribed to the detector response due to uncertainty in the calibration source activity and source-detector distance used in the detector response characterization. Systematic uncertainty associated with the source-detector orientation during the response characterization is expected to be small, but was not quantified here. Figure 2c shows the forward projection of the reconstructed PSL image (i.e., point-source and background) into data space overlaid on the measured data, summed across all detectors, as well as the deviance residual 42 between the measurement and fit. The fit represents the maximum likelihood solution given the point-source constraint implemented in PSL. The discrete separations observed in the deviance residuals are an artifact of low statistics and the discrete nature of the measured data (i.e., 0, 1 or 2 counts per bin). The attenuating effect of the steel containers is not included in the physics model and thus as the system surveys the back-side of the container stack, the forward-projected signal is expected to increase when approaching the source location. This causes two small unobserved peaks in the model (roughly 215 and 250 s into the measurement). The inclusion of attenuation in the reconstruction is part of an ongoing effort 43 but not included here. Finally, Fig. 2d shows an isometric rendering of the measured point cloud, colorized by the lidar return intensities as well as the PSL z-scores, demonstrating the rich 3D contextual detail that can be provided via lidar-based SLAM. While additional compute time is necessary for generating these visualizations, the reconstruction result was computed in approximately 27 s on the integrated GPU (Intel Iris Plus Graphics 650) on the on-board computer in MiniPRISM.
The Compton PSL results are similar to Fig. 2 and are included in Supplementary Fig. S1. They also include the true source location within 2 σ spatial bounds, and the best-fit voxel to true source location distance is 9 cm. The source activity estimate is 2.32 +0. 42 −0.40 (stat.) +0.22 −0.22 (sys.) mCi of 137 Cs (2σ uncertainty bounds). The estimate is consistent with the ground-truth activity as well as with the active coded mask PSL activity estimate. The total reconstruction time for the Compton PSL reconstruction was 7.7 s. The faster reconstruction time in this case is attributed to the list-mode implementation of the Compton PSL algorithm (see "Methods") and the smaller number of measured doubles events.
Distributed source imaging. Several hundred 1-µ Ci 22 Na point-sources, certified to 5% uncertainty by the manufacturer, were used to create a mock distributed source with known extent and activity. To ensure the imaging system was unable to distinguish the point-like nature of the sources, the individual sources were placed close together (5.08 cm pitch), the detector maintained at least a 50 cm standoff from the sources, and the voxelization scale of the reconstruction space was set at twice the source pitch.
The measurement took place in an active indoor laboratory environment wherein three distributed sources, separated by at least 5 m, were simultaneously staged. The distributions were constructed with the arrayed pointsources in distinct configurations to model radioactive spills and are shown in Fig. 3a. The first configuration contains a rectangular patch of sources ( 6 × 12 ) on top of a table as well as a square patch of sources ( 6 × 6 ) on the floor next to the table, with a total activity of 108 µ Ci 22 Na. The second configuration was a single smaller rectangular patch ( 6 × 4 ) on the ground next to a fume hood, with a total activity of 24 µ Ci 22 Na. The final configuration was an L-shaped patch comprised of three square patches ( 6 × 6 ), with two of the patches on a set of flat boxes and one patch up against a wall. The total activity of the final configuration was 108 µ Ci 22 Na.
The MiniPRISM system was hand-carried to survey each area of the laboratory in detail for 1-2 min to model an operation in which the operator had access to the real-time count-rate information. The total measurement time was 5.6 min and the generated map of the laboratory was ≈160 m 2 . The distances of closest approach to the center of the three source distributions were 0.53, 0.63, 0.54 m, in order of their description above and in Fig. 3a.
Similar to the point-source scenario, the results were generated as a single reconstruction including all measurement data. The reconstruction energy region-of-interest (ROI) was 511 ± 50 keV and a total of 39,183 singles and 4,379 doubles were recorded in the ROI over the entire measurement. Other reconstruction parameters include the following: 10 cm voxels, 125,754 occupied voxels (2.7% of total), 3,342 poses (10 Hz) and 10 MAP-EM iterations. An occupancy threshold of 10 points was used in order to reduce the number of occupied voxels and thus decrease the reconstruction time. A higher pose frequency of 10 Hz was used here in order to more accurately capture the higher-acceleration movements that can occur in hand-held operation compared to sUAS surveys.
The L 1/2 regularization coefficients for the active coded mask and Compton MAP-EM reconstructions were set to 2 × 10 −2 and 8 × 10 −3 , respectively. The coefficients were determined via a coarse grid search to better match the total activity reconstructions of the three sources compared to ground truth. Regularization strength is relative to the loss in the imaging objective function and, as a result, the optimal coefficients are expected to vary with statistics and imaging modality. The coefficients were determined specifically for this data and thus may not extend to other measurements, though more general and algorithmic optimizations of such factors are an active field of research 44  www.nature.com/scientificreports/ by summing all voxels in a 1.25 m radius sphere around the center of the true distributions. The radius was selected as the size at which all three source activity reconstructions plateaued ( Supplementary Fig. S2) and the ground-truth center point was used as a reference point to quantitatively compare the reconstruction with the ground-truth activity. The reconstructed activities of a.1 and a.3 (104.8 and 104.8 µCi) are in excellent agreement with ground truth of 108 µCi. The reconstruction of source a.2 is observed to be noisy and the activity estimate (11.4 µCi) underestimates the true source activity (24 µCi). This is believed to be because of a general lack of signal in this area, due to the lower overall source strength and the larger distance of closest approach in the survey. There is also activity attributed to the vertical surface of the cart that was slightly nearer to the survey trajectory. This incorrect location of activity nearer to the detector further skewed lower the source activity estimate. The active coded mask MAP-EM reconstruction is included in Supplementary Fig. S3. The reconstructions were similarly spatially accurate to the Compton results, though more diffuse, exemplifying the differences in imaging resolving power between active coded mask and Compton imaging with MiniPRISM. The reconstructed activities for sources a.1-a.3 were 118.2, 6.3 and 112.0 µ Ci of 22 Na, respectively, and the total reconstruction time was 140.8 s.
In contrast to the point-source results presented above, the distributed source reconstructions do not yet quantify uncertainty in the activity estimates due to spatial reconstruction error, however we ascribe a 9.4% uncertainty in the absolute scale of the response function used for reconstruction (see "Methods"). Future work could include exploring other approaches to provide uncertainty estimates in reconstructed activities, such as Markov Chain Monte Carlo and approximate highest posterior density credible regions in MAP estimates 45 . These approaches introduce additional computational challenges and thus further work is necessary to develop a near real-time uncertainty quantification estimate for distributed sources.

Discussion
The results presented above demonstrate that the quantitative SDF approach is capable of accurately and precisely localizing point-sources and imaging distributed sources in both the spatial and activity domains. The images can be used to quickly find point-sources or map distributed sources in wide-area or cluttered environments as well as provide detailed contextual information to the user in order to quickly plan the next set of actions. Furthermore, the methods can be executed on the timescale of seconds to minutes, allowing results to provide actionable feedback during the survey. Additionally, the quantified images can be used to compute dose-rate maps in areas of interest surrounding the source which is critical in contamination assessment or dose-minimizing activities. While this work focused on gamma-ray reconstruction, similar approaches can be implemented to enable quantitative SDF for neutron sources.
All SLAM-based approaches will be limited by the uncertainties in the real-time pose estimates. While on the order of cm in most cases, the uncertainties may become significant in near-field imaging applications or in cases where high angular acceleration is applied to the system (e.g., fast twisting in hand-carried operations). When coupled to systematic errors in the detector response, either from insufficient physics models or incorrect simulation geometries, the reconstructions can introduce unwanted image artifacts. Additionally, some amount of solution degeneracy will also exist due to the finite angular resolution of the imaging system (e.g., 7-10 deg in MiniPRISM), the limited statistics in the data, and the potentially large number of degrees of freedom (i.e., voxels) in the model. These effects can make the reconstruction susceptible to over-fitting and thus care must be taken to choose appropriate types and strengths of regularization. Operator understanding, training, and realtime feedback can also help to mitigate over-fitting by ensuring sufficient spatial sampling and pose variation in important regions and visually inspecting that repeated passes do not substantially alter the solution.
The 511 and 662 keV energies reconstructed here represent a unique region in the energy spectrum where both active coded mask and Compton imaging have utility. Below a few hundred keV, the Compton imaging modality degrades as the number of measured doubles events diminishes, the separation between two-site interactions decreases (increasing the uncertainty in the scattering axis), and the effect of Doppler broadening (increasing the uncertainty in the scattering angle) becomes more pronounced 46 , in which case only active coded mask should be used. Conversely, above several hundred keV the number of photo-electric absorptiononly events in a single 1-cm 3 CZT detector decreases significantly and Compton events become more likely, reducing the efficacy of the masking material, in which case only Compton imaging should be used. Here we simply demonstrate the capability to use both modalities separately, but work has been done with mixed-event reconstruction 47 in the energy regime where both modalities are still applicable and could be applied to SDF in future work.
Finally, the binary voxel-occupancy model (occupied and free) currently implemented is limited to reconstructing gamma-ray emissions to surfaces observed by the lidar. A tri-state occupancy model to classify unknown space, in addition to free and occupied space, would enable the ability to reconstruct within enclosed volumes and other areas not directly seen by the lidar. This capability will complement the work towards inclusion of attenuation in the physics model to account for both observed and possibly unknown attenuating material in the scene (e.g., material inside cargo containers in addition to the steel enclosure), enabling tomographic reconstruction.

Methods
This section describes the MiniPRISM free-moving gamma-ray imaging system, detector response characterization, image reconstruction algorithms, and the uncertainty quantification approach used in this work.
Detector system. The MiniPRISM system is a free-moving multi-modal omnidirectional gamma-ray imager, consisting of a partially populated 3D grid (6×6× 4) of 1-cm 3 coplanar-grid CZT detectors 35 www.nature.com/scientificreports/ array of independently readout detector modules facilitates both the active coded mask and Compton imaging modalities, enabling omni-directional imaging of sources across a wide range of energies relevant to nuclear security, safety and safeguards (tens of keV to several MeV). The detector system is paired with the Localization and Mapping Platform (LAMP) 18 which includes the auxiliary sensor package (visual camera, lidar, and IMU) and single-board computer for real-time contextual and gamma-ray data processing. The entire system is powered with a single 98 W-hr Li-ion battery, typically offering about one hour of operation on a single charge. Battery life may fluctuate depending on the computational demand from image reconstruction and SLAM. The system is platform agnostic, meaning it is capable of being hand-carried or mounted to a vehicle such as an automobile, ground robot, or sUAS.
Detector response characterization. Quantitative gamma-ray emission reconstructions require characterization of the angular response of the detector system. Here we characterize the MiniPRISM system response through experimentally benchmarked simulations. The simulations were performed using the Monte Carlo particle-transport code Geant4 31-33 and enabled parametric study of the angular response functions over the relevant emission energy domain. The simulations track the interaction positions and energy depositions of source photons and secondary particles in the detector array. Energy depositions are blurred according to a measured energy-dependent resolution function. Particles are emitted as a uniform cone beam pointed at the detector, with an apex at the source position and an opening angle defined such that the cone fully encompasses the detector system. Events of interest are tallied and used to compute a response η in units of effective area. For a particular source simulation (i.e., a single direction and energy), the effective area is computed using where N s is the number of simulated particles, r s is the source-detector distance, θ is the cone opening angle, and N d is the number of tallied events in the detector. For singles imaging, N d includes only the events that deposit energy within the specified energy ROI in a single detector module. Therefore, in addition to energy and direction, η is computed for each individual detector in the array. It is the unique modulation pattern of η across the detectors (for a given source direction and energy) that can be leveraged in reconstructing the source. In Compton imaging, two-interaction event sequences (i.e., doubles) with a total energy deposition within the energy ROI are tallied across the entire detector array. The result is an energy-and direction-dependent η for the whole array. Statistical uncertainties are included in the simulated η to account for Poisson counting statistics (i.e., δN d = √ N d ). Additional systematic uncertainties will inevitably exist for the simulated results, arising from model inaccuracies such as geometry, material composition and gamma-ray detection physics (e.g., charge transport/collection in CZT). For example, a small dead layer beneath the anode is a well-documented phenomenon 48 in coplanar-grid CZT with a size roughly the same as the grid pitch. Here we include a 1.25 mm dead layer in simulation to accommodate the 1.05 mm inner grid pitch as well as the larger grid lines near the edges 49 . We also use a 225 µ m surface cut on the remaining five sides of the cubes to account for reduced sensitivity due to electric field non-uniformities from the surface and from other nearby detector modules/ electronics. The value of this cut was set empirically to better match the measured response. Additional systematic uncertainties are not explicitly accounted for in the calculation of η , though are discussed below.
Benchmark measurements were performed with a 4π detector scanning system designed and built at LBNL. The scanning system consists of a rotatable platform, on which the detector system is placed, and a rotatable arm, on which the radiation source is attached (Supplementary Fig. S4). The combination of the platform rotation around the azimuth and the arm rotation along the polar direction facilitates a source position anywhere in the 4 π detector coordinate system. The source holder arm is extendable, allowing for the source to be placed at a varying distances away from the system ranging from 44-143 cm.
Several independent scans were performed with relatively strong laboratory-scale sources (hundreds of µ Ci), with gamma-ray energies ranging from 60 to 1332 keV. Source directions were nominally defined with the Hierarchical Equal Area isoLatitude Pixelation of a sphere (HEALPix) library 50 50 for a visual representation. The source-detector distance was set at 110 cm and considered to be in the far-field (>10 times the roughly 7 cm spatial extent of the detector array). The source remained static at each position for 5 min, resulting in a total individual scan time of 1 hr. In contrast to equation,(1), the sources were assumed to emit isotropically and thus the measured effective area was computed using where B is the branching ratio of the decay of interest, t is the dwell time, A is the calibrated activity (decay rate) marked on the laboratory source, is the decay constant and T is the time elapsed since the source activity calibration. Here the dwell time t is much less than the half-life ( τ 1/2 = log(2) · −1 ) of the measured sources. A background subtraction was performed with a linear fit between the average counts in neighboring windows on both sides of the spectral region-of-interest. Neighboring window widths varied across the measured source www.nature.com/scientificreports/ energies, as some sources (e.g., 133 Ba) have closely-separated gamma-ray emission lines and care was taken to not include source photons in the background estimation. The largest sources of uncertainty in the calculation of η arise from uncertainty in A, r s and N d . The laboratory sources were calibrated by the manufacturer and reported to have a 3% total uncertainty (99% confidence). The MiniPRISM system was placed on the scanning system and aligned by eye to be in the center of the rotating platform and in line with the height of the rotating arm. Given the potential for misalignment in two joints of the system, a conservative uncertainty of r s = 110.0 ± 2.5 cm was assumed. Poisson counting statistics are included for the uncertainty on N d as √ N d . Finally, the rotating arm is known to exhibit varying degrees of sag along the polar angle due to the weight and extension of the arm away from the pivot point. This effect is pronounced at source directions that fall in line with one of the axes of symmetry of the detector array, where outer detectors should perfectly occlude inner detectors. With a small amount of source-arm sag, additional streaming paths become available to many of the previously occluded detectors and the detection sensitivity can increase significantly, especially at lower energies. To mitigate the discrepancy between simulation and measurement, the simulations were performed with a 5-deg sag in the polar angle for source directions around the equatorial plane (HEALPix indices 4-7 for N side = 1). Figure 4 shows the singles and doubles effective area characterization results at 662 keV. The top pane shows the per-detector agreement of the single effective area of measurement (red) and simulation (blue) across the 12 source directions, where uncertainties are represented by the extent of each rectangular marker. The general trends and amplitudes are seen to qualitatively agree quite well. The middle pane illustrates the deviations with a more quantitative approach, displaying the differences of the per-detector effective areas between measurement and simulation (in units of Gaussian z-scores) as standard box-plots (light green) across the 12 source directions. Violin plots (light blue) are overlaid to provide a sense of the distribution of the data points along the vertical dimension.
The majority of simulations of individual detector efficiencies agree with the observed data to within the 1-2 σ , and the overall simulated detector system singles efficiency is within 1 σ for all measured angles, though some variance is observed for several source directions. First, directions 0-3 are positioned at a polar angle MiniPRISM far-field 662-keV gamma-ray detector response characterization (singles and doubles) at 12 source directions around 4 π defined by the HEALPix discretization scheme with N side = 1 . Measurements were performed with the scanning system described in the text, using a 182.5 µ Ci 137 Cs at 110 cm from the center of the detector array and 5-min acquisitions at each source direction. Simulations were performed using Geant4. The top pane shows the measured (red) and simulated (blue) full-energy singles effective area (mm 2 ) of individual detector elements, with shading denoting uncertainties (discussed in the text), across the 12 source directions. The middle pane shows a box-and-whisker plot of the per-detector effective area differences between measurement and simulation in the top plot, expressed in units of Gaussian z-scores ( σ ). The boxes and whiskers describe the interquartile range (IQR) and 1.5×ICR, respectively, and points indicate outliers outside the extent of the whiskers. A red line is drawn in each box to denote the median of the data. Finally, a violin plot is overlaid to provide a sense of the distribution (via a 1D kernel density estimator) of the data across the vertical dimension. The bottom pane shows the doubles effective area differences between measurement and simulation of the entire detector array, in units of Gaussian z-scores. www.nature.com/scientificreports/ of 48 • , corresponding to a source above the system where the electronics and battery enclosure occlude the detector array. The Geant4 model of the enclosure includes the battery, single board computer and analog-todigital converter (ADC) boards, though does not have significant detail to include other potential occluding and scattering media such as electronic components, fasteners and cabling. Therefore, it is expected that some level discrepancy exists for these source directions. The remaining source directions agree to 1-2 σ with much smaller variance, except for position 6 in which the lidar and camera are positioned directly in front of the detector array. Specific details of the lidar and camera internals and material composition are unavailable in manufacturer documentation and thus the Geant4 models are only approximate, leading to a larger variance in the effective area comparison at this position. Due to this finding, measurements with the MiniPRISM system are designed to reduce the time at which the lidar and camera package are placed between the suspected source location and the detector array. While not explicitly modeled or characterized here, we assume a human operator may have a similar effect when hand-carrying the system, and thus care is taken to also reduce the time at which the operator is between the detector and source. We acknowledge that this may not be possible in practice where potential source locations are not known a priori, and thus additional work is needed to refine the response in these directions. Finally, the bottom panel of Fig. 4 shows the differences between the measured and simulated doubles effective area for the entire detector array in units of Gaussian z-scores. Agreement is made for all source directions to ≤ 2σ.
In this work, quantitative SDF is first demonstrated for sources with gamma-ray energies of, or near, 662 keV. Therefore, the characterization of detector response at energies other than 662 keV will not be presented here. Following characterization, higher-resolution simulations ( N side = 16 or 3072 source positions) were performed at source emission energies relevant to the considered sources in this work (i.e., 511 and 662 keV). Example simulated responses are shown in Fig. 5 for singles (one detector) and doubles at 662 keV. Distinct anisotropy is observed due to occlusions from the auxiliary sensors, electronics and battery as well as self-occlusions from other detector modules in the array.
Imaging modalities. In addition to the gamma-ray source localization and mapping capabilities provided by signal modulation due to the motion of the system near an area of interest (i.e., proximity mapping), this work reconstructs gamma-ray images with the active coded mask and Compton imaging modalities. The active coded mask modality relies on the energy-and direction-dependent modulation of incident gamma-ray flux from both passive materials (e.g., battery, contextual sensors) as well as active materials (i.e., other detector elements) in the system 51,52 . Active elements acting as both detector and occluding material decreases the overall weight of the system by removing the need for heavy Pb or W masks, increases the detection sensitivity and imaging fieldof-view, and can be used to discriminate full-energy and down-scattered absorption events 53 . Individual detector modules are arranged such that the spatial pattern of events in the array uniquely (ideally) determines the direction of the gamma-ray source. However, as the energy and penetrative power of the gamma-rays increase (above a few hundred keV), the uniqueness of the full-energy signal modulation degrades and the likelihood of scattering events increases.
The design of the active detector array facilitates Compton imaging of gamma rays that scatter in one detector and are absorbed in another. Higher-order scattering sequences are excluded in this work as the analysis is significantly more complex and the Compton events are dominated ( 90%) by doubles events. The positions and energy depositions of the two interactions define the axis and opening angle of a cone that includes the direction of the incident gamma ray. The overlap of multiple cones in 3D reveals the source location.
Compton kinematics dictate that the cone generated by the two-interaction event sequence has a vertex at the first interaction location. In some imaging systems, the event time resolution and distance between interactions are such that time-of-flight (TOF) methods can be used to determine sequence ordering ( a → b or b → a) 54  www.nature.com/scientificreports/ For MiniPRISM, the distance between interactions is on the order of centimeters and the timing resolution is ∼2.4 µ s, therefore TOF cannot be used. If E ≤ 256 keV, kinematics assert that the first interaction will always deposit less energy than the second. However, the energy regime in which Compton imaging is primarily used here is > 256 keV and thus the sequence ordering can be ambiguous. Previous work has taken a single cone approach, assuming the higher energy deposition event occurred first 55,56 . To maximize sensitivity with limited statistics, our approach is to always project both possible cones with a simplified weighting for each cone based on the Klein-Nishina scattering and photoelectric absorption cross-sections. The weighting is similar to the event selection criterion used in 57 , though in that case still only one cone was projected.
In both active coded mask and Compton imaging, energy selection of photo-peak events is employed to image gamma rays from unique nuclear decay lines of specific radioisotopes of interest and to isolate source events from background. Furthermore, as our reconstruction approach is performed in a globally consistent adaptive image space, we can perform data cuts in both space and/or time (e.g., selecting a portion of the trajectory and associated radiation data near a suspected source location) to further isolate source from background. Scene data fusion. In the two source measurement scenarios described in the "Results" section and the Supplementary Information, the data collection to image reconstruction workflow is as follows. As the system is moved freely through an environment, the contextual sensors and SLAM solution provide real-time (up to 10 Hz) globally consistent estimates of the system trajectory (as a time series of position and orientation) and the surrounding 3D environment (as a point cloud). At a given point in time, the trajectory extents (with a configurable padding) are used to define the bounds on which the point cloud is voxelized at some set resolution (e.g., 20 cm). The generated voxel space retains the concept of occupancy -i.e., whether a voxel contains a point (or a configurable minimum number of points) from the point cloud, and the occupancy grid defines the reconstruction space.
The time correlation of radiation data and the system trajectory can be performed in two ways. As the radiation data is collected in list-mode (i.e., event-by-event), the system trajectory can be interpolated to the radiation event times, providing system pose estimates for each gamma-ray event or Compton sequence in the detector, facilitating list-mode type reconstruction algorithms. The radiation data can also be binned according to the system trajectory times in order to run bin-mode type reconstructions, though there is inevitable information loss with this approach as quick translations or rotations (e.g., twisting a hand-carried system) may not be captured given the frequency at which the system pose is computed. The nature of the data or the algorithm itself can sometimes dictate the type of reconstruction (i.e., bin-mode or list-mode). In cases where either option is available, there are computational advantages to choosing one over the other -e.g., bin-mode reconstruction times will scale as P × D where P and D are the number of poses and detectors, respectively, and list-mode will scale as the number of detected events, N. Apart from high count-rate environments, typically N < P × D and thus a list-mode implementation is favorable. In this work, the active coded mask PSL algorithm is implemented in a bin-mode reconstruction (due to current implementation constraints) while all other reconstruction algorithms (Compton PSL, active coded mask MAP-EM and Compton MAP-EM) are performed in list-mode. The pose-correlated radiation interaction data and 3D occupancy grid are then sent to the reconstruction algorithm to produce an image (i.e., gamma-ray emission estimates on the 3D occupancy grid).
Image reconstruction. The objective of the imaging problem as defined in this context is, given some measured data y , to form a non-negative image estimate consisting of J voxels ŵ ∈ R J + , with units of gamma-ray emissions per unit volume per second, that satisfies the following where b is a constant background rate, R(w) is a regularization applied to the image, β is a parameter to control the strength of the regularization, and ℓ is the negative logarithm of the Poisson likelihood where ⊙ denotes element-wise multiplication, Ŵ is the Gamma function, and ȳ is the mean estimate from our linear model ȳ = K ·ŵ +bt , where K is the system matrix that maps image space to data space and t are the measurement dwell times.
Here y is acquired in list-mode and consists of N events. In the singles reconstruction paradigm (i.e., proximity and active coded mask), each y n is a single interaction event in detector d at pose p with energy deposition E. The pose of the system represents the position x and orientation of the system in a global coordinate frame at time t (provided by SLAM). In the Compton reconstruction paradigm, each y n represents an event-pair (i.e., a Compton sequence) of two time-correlated interactions in detector pair d a and d b with energy depositions of E a and E b , respectively, at pose p with total energy deposition E = E a + E b , which is considered within the fullenergy window of the energy of interest. The two interactions are assumed to be a Compton scatter followed by a photoelectric absorption, though often we are unable to determine the sequence ordering ( a → b or b → a).
The {n, j} element of the list-mode system matrix K ∈ R N×J + effectively represents the probability of a gammaray emitted from voxel j producing event y n . In the case of singles imaging, this can be written as where η(r nj , � n , E) is the effective area of the detector at event n as seen from a source in voxel j with gamma-ray energy E, r nj is the vector connecting voxel j to the position of the detector at event n, n is the orientation of the detector at event n, ξ nj is the transmission probability of a source photon travelling along r nj , and ρ(r nj , z) is a regularized distance metric given by where z is a nominal size of the detector system included here to mitigate unwanted near-field amplification and | · | is the Euclidean norm. In the case of MiniPRISM, we set z equal to 10 cm. The transmission probability is always assumed to be unity here because of the additional computational burden of its calculation (i.e., raycasting) and we have yet to employ real-time methods to estimate the material composition of intervening occupied voxels, though this is an active area of research 43 .
In the Compton imaging case, K nj represents the intersection of both Compton cones associated with the y n event-pair with voxel j, where each cone is weighted by the Klein-Nishina and subsequent photoelectric crosssections for the sequence as well as the Compton effective area of event-pair n from a source in voxel j with total energy E = E a + E b .
To solve equation,(3), a form of the iterative Expectation Maximization (EM) approach is used. If no regularization is applied, the objective reduces to Maximum Likelihood (ML), in which traditional list-mode ML-EM 58 update equations can be applied. With regularization, we solve for the Maximum A Posteriori (MAP) estimate using list-mode MAP-EM updates where ⊘ denotes element-wise division, the sensitivity ζ ∈ R J + is the integrated response over all possible detected events (including those not measured in y ), and D is a vector of partial derivatives 59 . While fit in all bin-mode type reconstructions, background is currently not solved for in the list-mode formulation. The reconstructed weights will therefore be biased higher to account for background events in y, though in practice (and in the results presented here) we have observed this effect to be small. The inclusion of background fitting in list-mode reconstructions is left for future work.
Regularization is useful in mitigating over-fitting in the case of limited statistics and potentially large number of image voxels. We have found good performance with sparsity regularization to force low intensity voxels to zero, effectively removing the voxels from the reconstruction. Here we implement the L 1/2 regularization of the form used in Qian et al. 30 .
Finally, the size of the system matrix will grow as the dimensionality of the problem increases and, in some cases, can become too large to store in computer memory, necessitating an on-the-fly computation of system matrix elements when needed. However, parallelization of the forward (backward) projection operation is possible due to the independence of the voxels (events). To maximize the gain of this parallelization, the reconstructions are executed on a GPU. The implementation mitigates the memory burden as well as significantly increases computational performance, facilitating near real-time, order of seconds, reconstructions for most scenarios relevant in this application space.
Point-source localization. Gamma-ray source localization and mapping using free-moving imagers is susceptible to bias and over-fitting due to variability in sensitivity across the space, limited angular resolution and the general ill-posedness of the problem. Moreover, resultant spatial reconstruction errors can introduce large intensity estimation errors due to the inverse square falloff of counts with source-detector distance. As a result, we developed a low-dimensional method called Point-Source Localization (PSL) 28,29 to improve spatial and intensity reconstruction in the point-source domain. In this case, the source is constrained to a single image voxel and likelihood ratio tests are used to first compute source detection against a background-only model. If a source is detected, additional likelihood ratio tests are used to compute spatial confidence intervals around the maximum likelihood source position. The observed Fisher information matrix (i.e., the matrix of second derivatives of the negative log-likelihood function with respect to the point-source activity and background) is used to compute an approximate confidence interval for the source activity in each image voxel. The union of the intervals inside a predefined spatial confidence interval (e.g., 2 σ ) is used to compute a conservative estimate of the confidence interval for the maximum likelihood source activity. Additional details regarding uncertainty quantification can be found in Bandstra et al. 43 .

Data availability
The data used to generate the results will be made available on the Open Science Framework and access to the reconstruction software can be coordinated through the LBNL Intellectual Property Office.