Surrogate distributed radiological sources III: quantitative distributed source reconstructions

As described in Part II [2], in August 2021 we conducted a week-long experimental campaign at Washington State University (WSU) using radiation detectors borne on unmanned aerial systems (UASs) to measure the gamma-ray signature of various surrogate distributed Cu-64 source patterns. The source patterns were constructed out of arrays of up to $100$ point sources spaced densely enough to “look like” a truly continuous distributed source of Cu-64 for a given detector and trajectory (according to the method developed in Part I, Section II), and range in complexity from a uniform square source to regions of higher and zero activity superimposed on a uniform baseline (see Part I, Fig. 3). The nominal point source activities were ${\sim}7$ mCi at 0800 PDT of each measurement day, though the sources decayed throughout the day with a $12.7$-hour half-life, and the count rate comparisons made in Part II Section V-A suggest that the true source activities were ${\sim}1.4\times$ higher than the planned nominal activities, which were subject to substantial calibration uncertainties. In this work, we consider this ${\sim}1.4\times$ scale factor to be part of the overall system calibration. The predominant decay signature of Cu-64 is the $511$ keV annihilation photon line, which can be detected either in ‘singles’ mode (i.e., through its photopeak) or in ‘doubles’ mode (i.e., through its Compton scattering between two detector segments). In this work, we restrict ourselves to singles imaging, but note that Compton imaging analyses remains of interest. The detector systems used are NG-LAMP ($4\times 1^{\prime\prime}\times 1^{\prime\prime}\times 2^{\prime\prime}$ CLLBC modules) [9] and MiniPRISM ($58\times 1\,\text{cm}\times 1\,\text{cm}\times 1\,\text{cm}$ CZT modules) [10], both of which are capable of omnidirectional gamma-ray imaging. Due to the imperfect energy resolution of real detectors, we define the energy region of interest (ROI) corresponding to the $511$ keV photopeak to be $[481,581]$ keV for NG-LAMP and $[481,541]$ keV for MiniPRISM—see Part II. Also as described in Part II, when necessary we exclude data from detector crystals that regularly exhibited poor spectral and/or timing performance throughout the campaign. As such, $2/4$ NG-LAMP crystals and $53/60$ MiniPRISM crystals were used for these analyses. Finally, we restrict our analysis to data collected during the constant-altitude UAS raster patterns, i.e., we exclude takeoff/landing and any other free-form trajectories performed.

To make detailed quantitative comparisons between the distributed sources designed in Part I and their SDF reconstructions, accurate models of the measured scene must be generated and then aligned with the designed coordinate frame to as good a precision as possible.

The NG-LAMP and MiniPRISM detector systems are each equipped with a lidar and inertial measurement unit (IMU) to enable lidar-based simultaneous localization and mapping (SLAM) [11, 12] in order to reconstruct the detector trajectory and a 3D digital point cloud model of the mapped area. To co-register all SLAM results to the idealized field coordinate frame used in Part I, we first identified a high-quality SLAM map of a square source run in which most of the individual source locations were resolvable. The field had a slope of ${\sim}0.5^{\circ}$ for drainage, so we then aligned the field normal to the $+z$ axis by fitting a small plane near the center of the square source. After the point cloud was leveled, we were able to identify $64$ of the $100$ source locations through manual inspection in the CloudCompare software, and imputed a constant $z$ source position based on the apparent ground surface within the point cloud. These lidar source position estimates were then aligned with their corresponding points in the ideal square source (of exactly $4$ m grid spacing) using a least-squares distance minimization between all $64$ corresponding point sets, and the resulting transformation was applied to the leveled point cloud. The mean distance between the designed and lidar source points after co-registration was $6.2$ cm (with the mean $x$- and $y$-component differences $\ll$ 1 mm), and the mean grid spacing between the lidar source points was $4.0047$ m, only $0.12\%$ different from the design value of $4$ m. Fig. 1 shows the lidar point cloud colorized by height and the $xy$ positions of the $64$ discernible square source locations after alignment.

The remaining point clouds were then co-registered to the now-aligned reference cloud using the Iterative Closest Point (ICP) algorithm in Open3D [13]. We enforced a minimum ICP fitness score of $0.98$ (chosen empirically), and manually provided improved initial guesses for the transform if the initial ICP result did not meet this threshold on the first pass. Fig. 2 shows the results of this co-registration procedure on $35$ point clouds from the campaign. Finally, after all measurements were co-registered, the volume for radiological source reconstruction was defined as the plane $x\in[0,125]\cap y\in[0,80]$ m at $z=0$, with a $1\times 1$ m pixel size.

As described in Part I Section II-A, we perform single-energy gamma-ray imaging by assuming a linear model

where ${\boldsymbol{w}}\in\mathbb{R}^{K}_{+}$ is the vector of gamma-ray intensities or “weights” at each point or voxel center indexed $k=1,\ldots,K$, ${\boldsymbol{\lambda}}\in\mathbb{R}^{I}_{+}$ is the vector of expected counts observed by a detector over a series of measurements indexed $i=1,\ldots,I$, $\boldsymbol{V}\in\mathbb{R}^{I\times K}_{+}$ is the system matrix that relates the weights to the expected counts, and the $+$ subscript denotes non-negativity. An unknown but constant-in-time background term, $b$, can also be included and solved for (see Ref. [3]), but we exclude it in this discussion for clarity and because our measurements are heavily source-dominated. The $(i,k)$ element of the system matrix is

$\displaystyle V_{ik}=\frac{\eta_{ik}t_{i}}{4\pi|\vec{{\mbox{$\leavevmode\resizebox{6.50403pt}{5.78172pt}{\includegraphics[trim=10.00002pt 0.0pt 140.00021pt 0.0pt,clip]{scriptr/ScriptR}}$}}}_{ik}|^{2}}\exp\left(-\mu_{\text{air}}|\vec{{\mbox{$\leavevmode\resizebox{6.50403pt}{5.78172pt}{\includegraphics[trim=10.00002pt 0.0pt 140.00021pt 0.0pt,clip]{scriptr/ScriptR}}$}}}_{ik}|\right),$

(2)

where $\eta_{ik}$ is the effective area to source point $k$ during measurement $i$ (e.g., Fig. 5(a) of Ref. [3]), $\vec{{\mbox{$\leavevmode\resizebox{6.50403pt}{5.78172pt}{\includegraphics[trim=10.00002pt 0.0pt 140.00021pt 0.0pt,clip]{scriptr/ScriptR}}$}}}_{ik}$ is the corresponding distance vector, and $t_{i}$ is the $i\textsuperscript{th}$ dwell time. A measurement (or series thereof) consists of a Poisson sample

which has a corresponding Poisson negative log-likelihood of

where $\odot$ denotes element-wise multiplication. Eq. 4 can also be used in defining the “unit deviance” or “deviance residual” between the data and the model, which reduces to a $\chi^{2}$ metric at high counts—see Part I, § II-D.

A maximum likelihood radiation image reconstruction algorithm then seeks the optimization

or, by incorporating prior information, a maximum a posteriori algorithm seeks

for some hyperparameter $\beta>0$ and regularization function $f({\boldsymbol{w}})$, which are often chosen to promote sparseness and/or smoothness in the reconstructed image $\hat{\boldsymbol{w}}$. Eq. 5 can be solved iteratively through the maximum likelihood expectation maximization (ML-EM) algorithm [14, 15]

where $\hat{{\boldsymbol{w}}}^{(0)}$ is typically initialized to a flat image and the sensitivity map

describes the expected number of counts per unit weight for each source element $k$ summed over the $I$ measurements. Moreover Eq. 7 can be extended to solve the regularized minimization in Eq. 6 (“maximum a posteriori expectation maximization” or MAP-EM)—see Ref. [3], Eqs. 7–8 for further details. In particular we consider the sparsity-promoting, non-convex $L_{1/2}$ prior [16],

and the smoothing and edge-preserving total variation (TV) prior [17, 18], which can be written in 2D as

Computationally, the integrals in Eq. 10 are replaced with discrete sums and the derivatives with finite differences. An additional parameter $\varepsilon$ is often introduced to stabilize calculations when taking the gradient of Eq. 10 [18, §II-C]. Reconstruction algorithms are implemented in the mfdf [19] package using radkit libraries [20], and run on a graphics processing unit (GPU) via PyOpenCL.

To evaluate the quality of a reconstructed radiation image ${\boldsymbol{\hat{w}}}$ compared to its corresponding ground truth image ${\boldsymbol{w}}$, we will typically consider three metrics: (1) the ratio of total reconstructed and true activities,

(2) the normalized root-mean-square error (NRMSE),

and (3) the structure coefficient from Eq. 10 of Ref. [21],

for some small stabilizing constant $\epsilon>0$. We note that the structure coefficient, $s$, is essentially the Pearson correlation coefficient between the image pixel intensities. The activity ratio, $R_{\text{tot}}$, is chosen to give an overall intensity accuracy (closer to $1$ is better), the NRMSE is chosen to measure the average intensity similarity over the image independent of pixel position and normalized by the true total (closer to $0$ is better), and the structure metric, $s$, is chosen to measure perceived image similarity (closer to $1$ is better). In Figs. 3, 4, 7–10, and 12–14, these metrics are displayed on each reconstructed activity image.

Because the on-field source distributions consisted of ${\leq}100$ individual point sources specifically configured to emulate continuous distributed sources, for the purposes of quantitative image comparison the “true” images ${\boldsymbol{w}}$ are generated by interpolating the point source arrays to $1$-m pitch continuous ground truth images with the same average activity concentrations.

While reconstructed images are computed over the entire $125\,\text{m}\times 80\,\text{m}$ field, much of that area has zero source present, and, as shown in Section III, correspondingly quite low reconstructed activity compared to the areas with source present. Thus, to reduce the extent to which image quality metrics are dominated by near-zero comparisons, we compute the image quality metrics over smaller subsets of the image space. In particular, we compute the NRMSE using only pixels with non-zero activity in the interpolated ground truth distributions, and compute the structure metric in the domain $x\in[40,120]\cap y\in[0,80]$ m.

Figs. 3 and 4 show radiation image reconstruction results from eight different source measurements with the NG-LAMP system at a nominal $6$ m above ground level (AGL), a raster pass spacing of $5.2$ m, and a speed of $2.6$ m/s. Binmode MAP-EM reconstructions with a time binning of $t_{i}=0.2$ s were performed in a 2D image plane of dimensions $125\,\text{m}\times 80\,\text{m}$ with $1\,\text{m}\times 1\,\text{m}$ pixels, and ran in $1$–$2$ s on a 2019 MacBook Pro using an AMD Radeon Pro 5600M 8 GB GPU. Here we use the sparsity-promoting $L_{1/2}{}$ regularizer with coefficient $\beta=3\times 10^{-3}$ and $30$ iterations, which are useful representative values chosen manually based on our prior reconstruction experience. We note however that we explore the impact of these regularization parameters, as well as many of the other measurement parameters, in the subsequent sections. We also note that in Figs. 3, 4, and many subsequent image reconstruction figures, the measured trajectories are overlain and colorized by gross count rates (the aforementioned breadcrumbing approach) with a cyan-to-magenta color scale as a qualitative, additional guide, without showing explicit quantitative colorbars. In general, though, the gross count rates in these figures are ${\sim}1000/$s far from the source and as high as ${\sim}15\,000/$s above high activity concentration areas and/or sources measured early in the day. Both trajectories and point source array distributions are shown when their inclusion is relevant and they can be overlain without compromising the clarity of the reconstructed image.

In general, the reconstructions accurately map the shape and, after calibration, the magnitude of the ground truth radiation distributions. Some overall differences in activity distribution are notable, however. The interior activities of the shapes are generally over-predicted compared to ground truth, while the edges are under-predicted. In addition, the corners and edges of the reconstructed images are rounded compared to the sharp features of the true images, falling off smoothly over a few meters rather than instantly over a $1$ m pixel. As a result, the zero-activity corridor is not very sharp in the $8$ m separation reconstruction, though it is better-resolved in the $12$ m separation reconstruction. Quantitatively, the ratios of total reconstructed activity $R_{\text{tot}}$ range from $0.859$ ($8$ m separation) to $0.919$ (plume), indicating high precision and a slight ${\sim}10$–$15\%$ bias in the absolute magnitude of the source reconstructions. The NRMSE ranges from $0.151$ (hot/coldspot) to $0.380$ ($12$ m separation), indicating fairly good average shape agreement despite the smoother nature of the reconstructed images. Finally, the structure coefficient ranges from $0.870$ ($12$ m separation) to $0.935$ ($10\times 10$ square), indicating a high level of perceived shape similarity.

For additional context, Fig. 5 compares the measured and MAP-EM-reconstructed counts vs. time from the square source reconstruction of Fig. 3. The deviance residuals and smoothness of the MAP-EM curve qualitatively indicate that the MAP-EM algorithm is not overfitting to noise, and thus that $30$ iterations is a suitable choice in this scenario. For an example energy spectrum, see Fig. 5 of Part II. Finally, Fig. 6 shows a 3D render of the measurement and reconstruction.

Fig. 7 shows the results from a series of measurements in which the NG-LAMP system was flown over the $12$ m separation source at altitudes of $5$, $6$, $8$, and $10$ m AGL in order to observe trends in the reconstruction vs. detector altitude. All three performance metrics degrade monotonically with increasing altitude. The structure coefficient in particular falls from $0.883$ to $0.796$ and the sum ratio decreases from $0.885$ to $0.813$ between $5$ m and $12$ m AGL.

These trends confirm that in the idealized case of an omnidirectional radiation detector flying over a flat plane with no intervening material, there is no imaging advantage in flying higher, as opposed to a traditional optical camera where the size of the viewed area increases with altitude. In real measurement scenarios, there are of course operational advantages to higher altitudes, in particular, the avoidance of obstacles that may attenuate the radiation signal and pose a collision risk for UASs.

Fig. 8 demonstrates the reproducibility of the method, where the NG-LAMP system was flown over the $10\times 10$ square source at $6$ m AGL in four runs between approximately $1100$ and $1430$ PDT on the same day. The three performance metrics are remarkably consistent (less than a few percent variation) run-to-run, despite small variations in the measured trajectory and the decay of the Cu-64 sources ($t_{1/2}=12.7$ hours), both of which are visible in the images.

Fig. 9 shows the results of an analysis of data where the MiniPRISM system was flown over the $12$ m-separation source at a nominal $8$ m AGL and the fidelities of both the data and detector response models were progressively coarsened. In particular, the data were re-sampled from the original $t_{i}=0.2$ s to $t_{i}=5$ s and $10$ s dwell times, and the detector response was degraded from its initial $4\pi$ anisotropic multi-crystal response to a single isotropic detector with the same total effective area $\eta$. The coarsening of the response to an isotropic detector has only a small negative effect on the reconstruction performance; this is not surprising as the $r^{2}$ modulation of the counts in Eq. 2 is much stronger than the $\eta$ modulation for these airborne surveys. For further exploration of the negative effect of response coarsening on reconstruction performance, we refer the reader to Ref. [22]. The coarsening of the data to $5$- or $10$-s intervals, conversely, slightly and then significantly degrades reconstruction performance. While the degradation to $t_{i}=5$ s is small according to the three performance metrics, the shape of the distribution noticeably degrades by eye, with the two rectangular activity bands smearing towards the center of the image. At $t_{i}=10$ s, the degradation is even more substantial, with prominent activity smearing (concentrated near the trajectory points) and corresponding changes in the NRMSE and structure metrics.

Similar to the coarse-graining study of Fig. 9, Fig. 10 shows reconstruction results for the $12$ m-separation source with MiniPRISM at $8$ m AGL as the UAS flight speed is increased. Although during measurements only a single speed of $2.6$ m/s was set, we emulate faster speeds by downsampling the listmode radiation data and compressing the radiation and trajectory timestamps by a constant factor. In particular, Fig. 10 shows effective speed increases of $1$, $3$, $6$, and $9\times$, resulting in emulated speeds of $2.6$, $7.8$, $15.6$, and $23.4$ m/s, and corresponding total measurement times of $714$, $238$, $119$, and $79$ s. Note that we maintain a constant readout interval of $t_{i}=0.5$ s, unlike the coarse-graining study, and thus the baseline $1\times$ reconstruction in Fig. 10 differs slightly from the baseline reconstruction in Fig. 9.

As the speed of the UAS raster increases, the reconstruction noise increases and the image quality decreases due to the reduction in photon statistics. Fig. 11 shows the image quality metrics at more speeds than Fig. 10. While the image quality does not markedly decrease at a $3\times$ speedup ($7.8$ m/s), the NRMSE rises from the baseline $0.354$ to $0.478$ at the $6\times$ speedup ($15.6$ m/s), suggesting a top speed of ${\lesssim}8$ m/s is suitable for this particular mapping scenario. Fig. 11 also suggests that the NRMSE and structure metrics are strongly anti-correlated.

Fig. 12 shows the results of a series of analyses of a single data collection where the NG-LAMP system was flown at $6$ m AGL above the $\mathsf{L}$-shape source, and the $L_{1/2}$ regularizer coefficient and number of iterations were varied between relatively low and high values—$10^{-4}$ vs. $10^{-2}$ and $10$ vs. $50$, respectively. In all four cases, the overall $\mathsf{L}$-shape is correctly observed, but at various levels of sharpness and thus agreement with the true distribution. The $(10^{-2},\,50)$ image is over-regularized, producing a distribution that is significantly smaller than the true source distribution, and has a correspondingly high NRMSE of $0.832$ and low structure coefficient of $0.628$. At the other end of the spectrum, the $(10^{-4},\,10)$ image has not fully converged, resulting in a modest NRMSE of $0.388$ but an improved structure coefficient of $0.870$. The intermediate case of the $(10^{-4},\,50)$ image performs the best of the four combinations, with an NRMSE of $0.327$, a structure coefficient of $0.903$, and a sum ratio of $0.947$. This hyperparameter selection even slightly outperforms the initial hyperparameters $(3\times 10^{-3},30)$ used for the $\mathsf{L}$-shape reconstruction in Fig. 3, which had an NRMSE of $0.357$, a structure coefficient of $0.896$, and a sum ratio of $0.868$. Therefore we note again that the question of rigorous optimal hyperparameter selection based on comparison to ground truth is left for future work.

Fig. 13 shows the results of a second regularizer study in which NG-LAMP was flown at $6$ m AGL over the hot/coldspot source, and the TV regularizer coefficient $\beta$ and number of iterations were varied from $10^{-7}$ to $3\times 10^{-7}$ and $100$ to $1000$, respectively, while using a fixed $\varepsilon=2.4\times 10^{5}$ $\gamma/$s as per Ref. [18, §II-C]. The upper bound of the regularizer coefficient was kept lower than the $L_{1/2}$ regularizer study, and the final two image iterations were averaged together, in order to avoid inducing strong “checkerboard” image artifacts that oscillate every iteration at $\beta=3\times 10^{-7}$. Such checkerboard artifacts appear to be a limitation of the TV regularizer at relatively high $\beta$ arising from the fact that minimizing finite differences in $x$ and $y$ smooths the horizontal and vertical bands of the image, but does not constrain diagonally-adjacent pixels [23]. We also note that for computational tractability, these reconstructions use data summed over all NG-LAMP detector crystals. Both images using $100$ iterations appear under-converged, with some apparent blurring between the three hot, cold, and baseline regions. After $1000$ iterations, the hot and coldspots become more distinct from the uniform baseline, with the $\beta=10^{-7}$ image outperforming the $\beta=3\times 10^{-7}$ image in terms of NRMSE and structure coefficient. In all four TV-regularized images, the summed activity ratio is in general closer than images using the $L_{1/2}$ regularizer by about $10\%$, though the reason for this is not readily apparent, and similarly the NRMSE values are lower by up to $25\%$. Further discussion of the TV vs $L_{1/2}$ regularizer performance is given in Section IV.

Fig. 14 shows the results of another series of analyses from a single data collection in which NG-LAMP was flown at $6$ m AGL over the $\mathsf{L}$-shape source using a tighter-than-normal line-spacing of $4.1$ m (instead of $5.2$ m) to facilitate testing trajectory cuts in post-processing. In particular, we cut the trajectory into individual raster lines by simple timestamp cuts and then perform reconstructions (both unregularized and with the standard $\beta_{L_{1/2}{}}=3\times 10^{-3}$) using only every $n\textsuperscript{th}$ raster line for $n=1,\ldots,5$.

The correct shape of the $\mathsf{L}$-shape source is visible for $n=1,\,2,\,3,$ but the shape severely degrades for $n=4,\,5$. The degradation is more severe in the $L_{1/2}$-regularized reconstructions due to the sparsity-promoting nature of the regularizer. These trends are borne out by the overall drop in structure coefficient from $0.902$ at $n=1$ to $0.519$ at $n=5$ (regularized) and from $0.896$ to $0.827$ (unregularized), and the corresponding increase in NRMSE from $0.342$ to $0.879$ (regularized) and $0.329$ to $0.404$ (unregularized). For reference, Fig. 14 also shows the sensitivity maps for each post-cut measurement. The worst-performing images are associated with sensitivity maps that contain noticeable visual banding between the raster lines, indicating that achieving uniform sensitivity over the source area improves quantitative reconstruction. A similar line of reasoning can also be used to argue that some non-trivial amount of measurement time should be spent flying over source-free areas in order to better capture the contrast between signal and background, though we do not take up this question quantitatively here.