Localizing hot spots in Poisson radiation data matrices: nonnegative tensor factorization and phase congruency

Detecting and delineating hot spots in data from radiation sensors is required in applications ranging from monitoring large geospatial areas to imaging small objects in close proximity. This paper describes a computational method for localizing potential hot spots in matrices of independent Poisson data where, in numerical terms, a hot spot is a cluster of locally higher sample mean values (higher Poisson intensity) embedded in lower sample mean values (lower background intensity). Two numerical algorithms are computed sequentially for a 3D array of 2D matrices of gross Poisson counts: (1) nonnegative tensor factorization of the 3D array to maximize a Poisson likelihood and (2) phase congruency in pertinent matrices. The indicators of potential hot spots are closed contours in phase congruency in these matrices. The method is illustrated for simulated Poisson radiation datasets, including visualization of the phase congruency contours. The method may be useful in other applications in which there are matrices of nonnegative counts, provided that a Poisson distribution fits the dataset.

• quantifying the relative positions of the clusters.
We discuss a numerical computational method for this objective. The method is a combination of commonplace algorithms (for example, the discrete Fourier transform) and relatively newer results in tensor math for numerical data.
The numerical output of this processing is convenient for visual display and for superimposing on other digitized data, such as registered visual images (cf. [11]) or geospatial maps of a region of interest. For example, ATAK [12][13][14] is an application for Android smartphones and tablets which uses GPS together with maps of an operational area to display real-time information to field personnel. GPS-registered localization of potential radiation hot spots can be transmitted from a centralized server to ATAK devices and superimposed on infrastructure maps for local situational awareness. Near real-time numerical computation is expected on fast, inexpensive Graphical Processing Units [15].
"Tensor factorization and phase congruency for 2D frames of Poisson data" section describes the Poisson assumptions about the digitized sensor data and introduces the two key numerical computations. The computations are discussed and illustrated by simulation in "Nonnegative tensor factorization: Poisson data" and "Phase congruency in a 2D grid" sections. "Conclusions" section is conclusions. The Appendix has succinct summaries of the two computations, with references that contain the algorithmic details.

Tensor factorization and phase congruency for 2D frames of Poisson data
Various systems output sensor data at discrete time stamps (cf. [2-4, 16, 17]) with different physical interpretations and different digital formats. This paper concentrates on digitized sensor data in the format of 2D matrices in which the data is total measured counts in a Poisson process [18]. An example of contemporary technology is a coded aperture system for which an output matrix of Poisson data is computed when the input data is Poisson (cf. [19][20][21][22][23]).
We will call the 2D data matrix a frame. Each frame is an X × Y grid of independent Poisson counts. When a sequence of frames is generated over time, each frame has a discrete time stamp or index marking its place in the sequence. We assume that the set-up for data-acquisition is the same for every frame in a sequence, in particular, the acquisition time interval for recording gross counts is the same and the points in the 2D grids remain registered with respect to the object or area being scanned.
Given n frames with individual time stamps t 1 < t 2 < · · · < t n , the 3D dataset D is the X × Y × n array in which the n frames are "stacked" in order 1, 2, ..., n. By assumption, background radiation is independent and identically distributed ( iid ) Poisson counts with fixed intensity B at each grid point in D . A potential hot spot in D is a cluster or subgrid H of iid sample values which are greater on average than the iid sample values at neighboring grid points. Unless stated otherwise, it is assumed that the location, shape, and Poisson intensity H of a hot spot H do not change in a sequence of frames. (Unshielding a hot spot and small shifts in position are discussed respectively in the short "Unshielding a source" and "Real or apparent motion of a hot spot" sections.) Two numerical computations constitute our two-step procedure for detection and delineation of potential hot spots in dataset D : 1. nonnegative tensor factorization of 3D array D [24,25]. 2. phase congruency in pertinent 2D arrays [26][27][28].
The tensor factorization is computed to maximize a Poisson likelihood conditioned on the sample dataset D . This yields a second 3D array M to augment D ; then phase congruency is computed for both D and M . The coordinates of closed contours in phase congruency define the position and approximate outline of potential hot spots in D. Clusters in scatter plots of phase congruency are additional information.
Phase congruency is a frequency-based analysis based on well-known methods in discrete math, including wavelets and discrete Fourier transforms. The tensor factorization is a relatively new development in tensor math. Computation with tensors is more complex than with conventional 2D matrices, but the tensor factorization can reveal numerical relationships in higher dimensions.

Notation and computed tensor M
Given the nonnegative integer array D , we compute a nonnegative real-number array M to maximize the conditional Poisson likelihood P(D|M) . The objective is localization of potential hot spots in D , to which end the computed array M augments the dataset D both for visual comparisons and for the computations in "Phase congruency in a 2D grid" section.
The array M is computed in tensor math. In computational linear algebra, an nth-order tensor is an n-dimensional array of real numbers [29][30][31]. Thus, the 64 × 64 × 20 dataset D in the following "Simulation" section is a 3rd-order tensor with indices 1 to 64 in its tensor mode-1 and mode-2 (the first and second dimensions in the 3D array) and index 1 to 20 in its tensor mode-3 (the third dimension).
Computations in tensor math that minimize certain error functions have become important tools in applications involving real-number arrays of data (cf. [24,25,29,32,33] and their references). "Nonnegative tensor factorization of Poisson data" section in Appendix is a succinct outline of the computation of the tensor M to maximize P(D|M) by solving an equivalent minimization problem [24,25].
The computed tensor M has the same dimensions as the data tensor D . For short notation, let index i denote a 3D-location in tensor D. Let The sample variance of M in each of its modes is smaller than the corresponding sample variance of D: This tends to smooth the values in subgrids of similar values in M while sharpening the boundaries between adjacent subgrids of higher and lower values. Small subgrids or point sources with insufficient sample count in D are not guaranteed to emerge strongly in M , especially if dominated by nearby subgrids with higher total counts. Overall, however, recurrent subgrids are accentuated more in M than in D , and this is an advantage in the computation of phase congruency in "Phase congruency in a 2D grid" section. The optimal tensor M is computed in a factorized form for which the user must specify the number of components R [24,25,34]. We use R = 10 , a number found empirically for datasets in this paper such that changes in the numerical values computed in M are relatively insignificant for R > 10.

Simulation
A dataset D is simulated as follows: • Each frame is a 64 × 64 grid. Background radiation is simulated for each frame separately as The sample mean is µ D = µ M = 0.1072 and the sample variances are

Phase congruency in a 2D grid
The method of phase congruency Both the data tensor D and the optimal Poisson-likelihood tensor M in "Nonnegative tensor factorization: Poisson data" section are discrete grids in which a potential hot spot is a subgrid of higher mean value embedded in a neighborhood of lower mean value. A method based on Fourier frequencies is used for automatic (unsupervised, noninteractive) detection and delineation of these subgrids. In practice, an unsupervised method should accommodate different ranges of numerical values in grids and diverse subgrid shapes. Phase congruency in a matrix of real numbers does this by computing agreement or congruency locally in a 2D grid in the frequency domain [26][27][28].
In the Fourier representation of signals, including multidimensional arrays of numerical data, the unique role of phase in locating "events" such as edges or points has long been recognized (cf. [35]). In conventional image processing, phase congruency assigns   an invariant measure of significance to localized edges, lines, and corners [26][27][28][36][37][38]. For the numerical data in D and M , it reveals boundaries separating subgrids of higher average values from lower average values, thereby delineating clusters with higher sample means embedded in local neighborhoods of lower sample means. The coordinates of closed contours in 2D phase congruency define the location and outline of a potential hot spot cluster.
Cluster boundaries in 2D grids are step discontinuities characterized by coherence in phase of Fourier frequency components at several scales and orientations. Kovesi's refined algorithm [26] with noise compensation uses wavelets for local frequency information at a fixed number of scales and filters at a fixed number of orientations. Given an X × Y matrix, the computation returns an X × Y array of values in the range 0 to 1 where 0 indicates no significant "event" and 1 indicates high significance. Phase congruency is a dimensionless, normalized measure, and the information in phase congruency covariance matrices is conveniently displayed as contour and scatter plots.
In practice, the user must assign parameters such as the number of wavelet scales and the number of orientations. "Phase congruency in 2D grids" section in Appendix lists the values recommended in the literature [26] and used in this paper. Figure 3 shows phase congruency computed for the data D in Fig. 1 with the two hot spots H 1 and H 2 . Figure 3b is a contour plot for phase congruency in the 64 × 64 mode-3 sample mean values, the sample values averaged over the 20 frames. The 2D grid Clusters for H 1 and H 2 in the scatter plots (a), (c) are irregular and incomplete but align with the hot spot contours in (b); however, clusters are also seen outside the true hot spots. These extra clusters are created by random events in the background Poisson process, specifically, by random "Poisson clumps" [39] of grid points with higher sample values than the neighboring values. These events occur randomly from frame to frame, but when averaged over all 20 frames, most counts are too small to emerge from background as phase congruency contours in (b). Figure 4 shows phase congruency in the computed tensor M in the same layout as Fig. 3. The two closed contours in Fig. 4b correlate strongly with Fig. 3b; however, the  Fig. 1. a, c Scatter plots for the 20 frames along mode-3. b 3-level contour plot scatter plots in Fig. 4a, c have fewer random clusters and better delineation of the two hot spot clusters aligned with the contours in (b). The gaps in plots for H 2 in frames 1 to 7 in Fig. 4a and frames 3 to 6 in Fig. 4c are due to random samples in H 2 with comparatively low values in those projections at those particular frames. Figure 5 shows phase congruency in the data D in Fig. 2 with the three hot spots H 1 , H 2 , H 3 . Figure 6 has the same layout for the computed tensor M .

Unshielding a source
If a hot spot is heavily shielded by other material initially but the shielding is removed at some time stamp in a sequence of frames, clustering in mode-3 in tensors D and M facilitates detecting its unshielding. (Evolution might also occur, as in a sequence of frames of the uptake of a radiotracer by an organism [17]).
The data D in Fig. 7 has H 2 shielded in frames 1 to 9, then unshielded in frames 10 to 20. Figure 8

Real or apparent motion of a hot spot
There may be real or apparent movement of a hot spot in a sequence of frames. Repositioning a sensor array creates an apparent shift of a stationary source and misaligns with previous frames, but if the real or apparent motion is minor relative to the generation of frames, then the hot spot shift is small in consecutive frames in D. Figures 9 and 10 show a small shift in H 2 at frame 10. The dotted-line arrow in the contour plot (b) indicates the direction of the shift. The scatter plots (a), (c) show the relatively small shift in the H 2 cluster beginning at mode-3 index 10 (frame 10), with fewer random clusters and better delineation in Fig. 10 than Fig. 9.

Conclusions
This paper presents a two-step, numerical computation for automatically localizing potential hot spots in matrices of gross Poisson counts: 1. Given a 3D dataset D of 2D frames of Poisson count data, a 3D tensor factorization M is computed to maximize the conditional Poisson likelihood P(D|M) . The maximization of this likelihood is achieved by minimizing a Kullbach-Leibler divergence function [25,40], specifically, M is computed in tensor math to minimize the function i m i − x i log m i for index i over the nonnegative integers x i in D and the corresponding nonnegative reals m i in M. 2. Phase congruency is computed for the two grids D and M projected in their three tensor modes. Phase congruency provides an invariant, normalized, numerical measure of "events" in D and in M . Kovesi's method [26] incorporates compensation for noise, a weighting for the spread of local frequencies, and filtering attuned to step discontinuities in diverse orientations in the grids.
The computed tensor M augments the data tensor D for analysis of phase congruency. There is strong correlation between D and M in the contours of phase congruency in the sample mean values over a sequence of frames, but there are greater differences in the contours of phase congruency in the sample means along the axis of frames (along the capture-time index). Random clusters may occur in a sequence of frames in D due to Poisson clumping of background counts, even at a low Poisson rate. The computation of M tends to suppress these random clusters and reveal the contours of recurring clusters. Areas of current research and development include: 1. The mathematical relationship between background Poisson probabilities and the level sets of phase congruency in a grid D . The goal is a probability-based decision rule for potential hot spots based on results in Poisson Scan Statistics (cf. [41][42][43]). This work includes a characterization of false positives (random clusters incorrectly called hot spots) and false negatives (failures to detect true hot spot clusters). Bayesian Spatial Scan Statistics [44] may yield lower error rates. 2. Fusion of data from multiple sources. Some applications have multiple radiation monitoring systems [9] or multi-modal sources of information [45]. In certain situations, tensor math facilitates the fusion of numerical data [46] into a composite situational map for a geospatial area.
The application in this paper is matrices of Poisson radiation data; however, the method is potentially useful in other applications involving matrices of nonnegative integer counts that are described by a Poisson probability law. We point out that if the average counts are small and the occurrences of high counts are rare events in a dataset, then a Poisson approximation might be justified by the so-called "Poisson law of small numbers" [47].

Nonnegative tensor factorization of Poisson data
Background references in tensor computation with definitions, basic tensor math, and aspects of numerical algorithms include [29][30][31][32][33]. This section summarizes the nonnegative tensor factorization of a data tensor D of independent Poisson samples. In contemporary numerical math, an nth-order tensor is an n-dimensional array of real numbers. In this terminology, a 1D vector is an order-1 tensor and an M 1 × M 2 matrix is an order-2 tensor. If tensor D is an M 1 × M 2 × · · · × M n array, then its mode-1 index is the integer range 1 : M 1 , its mode-2 index is 1 : M 2 , ..., its mode-n index is 1 : M n .
Classic matrix computation is a starting point for developments in tensor computation [30]. In particular, many aspects of Singular Value Decomposition (SVD) of matrices generalize to Higher-Order SVD (HOSVD) of nth-order tensors [29,33]. Further results in tensor math concern factorizations for least-squares error (LSE); specifically, given tensor D , the tensor D is computed in a designated factorized form to minimize the norm ||D −D|| F [32,33]. (The Frobenius norm ||Y|| F of tensor Y is the square root of the sum i y 2 i over all the values y i in Y [30]). Implicit in using LSE for a rank-1 tensor factorization of D is an assumption of Gaussian variation in the data [24]. An alternative is to compute a Poisson model M for tensors of nonnegative counts which might be sparse [24,25]. This includes nonnegative arrays of radiation data D which is Poisson in distribution.
Given order-n tensor D of nonnegative integer counts, the tensor M is computed to maximize the conditional Poisson likelihood P(D|M) . M is a tensor of nonnegative real numbers the same size as D . Letting index i denote a 3D-location in D and x i [or The iterative computation developed in [24,25] computes the optimal tensor M in a rank-1 factorized form. For an nth-order tensor D and a user-specified number of components R, the method computes scalors α r for 1 ≤ r ≤ R and matrices A (i) of respective sizes M i × R for 1 ≤ i ≤ n ; then the tensor M is where a (i) r is the rth column of matrix A (i) and • is outer product. The outer product of n matrices creates an nth-order tensor.
Convergence of the iteration to a tensor M maximizing P(D|M) has been proved for mild conditions on the nonzero values in D (roughly, both the density of nonzero values and their spread with respect to the size of the component matrices A (i) must be adequate) [24]. The implementation in MATLAB Tensor Toolbox version 2.6 [34] is a version of the alternating Poisson regression in [25]. We used default values of the control parameters (see [34]) except that the stopping tolerance was 0.5e−03. Random number generation was reset to its default seed at the start of each run so that the iteration began with the same "randomized" initial tensor. As noted in "Nonnegative tensor factorization: Poisson data" section, the number of components R = 10 for these datasets is an empirical value where, for R > 10 , the computed tensor M has relatively insignificant changes in numerical values.
Tensor factorization reveals multilinear relationships in a multidimensional numerical dataset. There are well-documented numerical methods to maximize conditional Poisson likelihood without tensor math (cf. [21-23, 38, 40]), but a mathematical model of the source-to-detector projection of data is required.

Phase congruency in 2D grids
Phase congruency in 2D finds boundaries between subgrids of higher mean values and adjacent neighborhoods of lower mean values, or vice versa. When the grid is sample mean values from a spatial Poisson process like dataset D or from an array of Poisson parameters like computed tensor M , the boundaries delineate clusters of samples iid at higher Poisson intensities, or potential hot spots, embedded in background samples iid at a lower Poisson intensity. See [26][27][28] for details of phase congruency in 2D using wavelets. This section is a high level description and includes the control parameter settings recommended in the literature [26]. The input matrix is an X × Y grid. For each grid location, and each scale and orientation, the initial energy is found by convolving the 2D array with even-odd quadrature filters, then the results are processed in these steps: 1. Randomness in phase due to noise is estimated. Noise-suppression parameter k = 2 indirectly sets a reference for noise in phase that is subtracted locally. 2. To insure that the spread of frequencies is adequate, a sigmoid weighting function penalizes too few in-phase frequencies. Sigmoid function parameter c = 0.5 is the cut-off below which deemphasis kicks in, and parameter g = 10 controls the sharpness of the function. 3. The adjusted energies are summed over all filter orientations and divided by the sum of the response amplitudes of the wavelets over all orientations and scales. This paper uses 4 wavelet scales and 6 filter orientations in 2D. To limit the spatial extent of the phase analysis in 2D, frequencies with wavelengths larger than min = 3 are suppressed.
The maximum and minimum moments of phase congruency covariance are X × Y matrices of phase information from all scales and orientations [28]. These matrices are computed in a suite of downloadable MATLAB programs [26] and their average is used in this paper. Phase congruency has been implemented for 3D grids [48] like dataset D ; however, the 3D computation and visualization are more complicated than 2D, and a preliminary evaluation did not show substantial changes in the outcomes for the 3D datasets in this paper.