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

Abstract

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.

Introduction

Applications of radiation monitoring systems range from the detection of unsecured radioactive sources in otherwise normal environments to situational awareness after a nuclear disaster. Fixed-position sensor systems are popular (cf. [1, 2]) and there is increasing use of mobile sensors—ground-based [3,4,5,6], on aerial platforms [6,7,8], or combinations of both [9, 10].

A two-stage radiation monitoring system has acquisition and digitization of sensor data in a first stage followed by data processing in a second stage (cf. [10]). This paper focuses on data processing where the digitized sensor data is matrices of Poisson counts. In many applications, an objective in the processing stage is detection and delineation of potential radiation hot spots which, in numerical terms, means

• discovering clusters of comparatively high radiation intensity in the digitized sensor data, and

• 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,3,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\times 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< \cdots < t_n$, the 3D dataset ${\mathbf{D }}$ is the $X\times Y\times n$ array in which the n frames are “stacked” in order 1, 2, ..., n. By assumption, background radiation is independent and identically distributed (${\mathbf{iid }}$) Poisson counts with fixed intensity $\lambda _B$ at each grid point in ${{\mathbf{D }}}$. A potential hot spot in ${\mathbf{D }}$ is a cluster or subgrid H of ${\mathbf{iid }}$ sample values which are greater on average than the ${\mathbf{iid }}$ sample values at neighboring grid points. Unless stated otherwise, it is assumed that the location, shape, and Poisson intensity $\lambda _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 ${\mathbf{D }}$:

• 1. nonnegative tensor factorization of 3D array ${\mathbf{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 ${\mathbf{D }}$. This yields a second 3D array $\mathbf{M }$ to augment ${\mathbf{D }}$; then phase congruency is computed for both ${\mathbf{D }}$ and $\mathbf{M }$. The coordinates of closed contours in phase congruency define the position and approximate outline of potential hot spots in ${\mathbf{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.

Nonnegative tensor factorization: Poisson data

Notation and computed tensor M

Given the nonnegative integer array ${\mathbf{D }}$, we compute a nonnegative real-number array $\mathbf{M }$ to maximize the conditional Poisson likelihood $P({\mathbf{D }}|\mathbf{M })$. The objective is localization of potential hot spots in ${\mathbf{D }}$, to which end the computed array $\mathbf{M }$ augments the dataset ${\mathbf{D }}$ both for visual comparisons and for the computations in “Phase congruency in a 2D grid” section.

The array $\mathbf{M }$ is computed in tensor math. In computational linear algebra, an $n{\text{th}}$-order tensor is an n-dimensional array of real numbers [29,30,31]. Thus, the $64\times 64 \times 20$ dataset ${\mathbf{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 $\mathbf{M }$ to maximize $P({\mathbf{D }}|\mathbf{M })$ by solving an equivalent minimization problem [24, 25].

The computed tensor $\mathbf{M }$ has the same dimensions as the data tensor ${\mathbf{D }}$. For short notation, let index i denote a 3D-location in tensor ${\mathbf{D }}.$ Let $x_i$ [or $m_i$] denote the value in ${\mathbf{D }}$ [or $\mathbf{M }$] at location i. A value $m_i$ in $\mathbf{M }$ is usually not a simple sample mean of the data ${\mathbf{D }}$ in any single dimension; however, the sum $\sum _i m_i$ in $\mathbf{M }$ equals the sum $\sum _i x_i$ in ${\mathbf{D }}$, and the sample mean value $\mu _{{\mathbf{D }}}$ of the entire dataset ${\mathbf{D }}$ is preserved in the projection of tensor $\mathbf{M }$ in each of its three modes:

$$\mu _{\mathbf{D }} = \mu _{\mathbf{M }_1} = \mu _{\mathbf{M }_2} = \mu _{\mathbf{M }_3}.$$

The sample variance of $\mathbf{M }$ in each of its modes is smaller than the corresponding sample variance of ${\mathbf{D }}$:

$$\sigma ^2_{\mathbf{M }_i} < \sigma ^2_{{\mathbf{D }}_i} \text{ for } i=1,2,3.$$

This tends to smooth the values in subgrids of similar values in $\mathbf{M }$ while sharpening the boundaries between adjacent subgrids of higher and lower values. Small subgrids or point sources with insufficient sample count in ${\mathbf{D }}$ are not guaranteed to emerge strongly in $\mathbf{M }$, especially if dominated by nearby subgrids with higher total counts. Overall, however, recurrent subgrids are accentuated more in $\mathbf{M }$ than in ${\mathbf{D }}$, and this is an advantage in the computation of phase congruency in “Phase congruency in a 2D grid” section.

The optimal tensor $\mathbf{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 $\mathbf{M }$ are relatively insignificant for $R > 10$.

Simulation

A dataset ${\mathbf{D }}$ is simulated as follows:

• Each frame is a $64\times 64$ grid. Background radiation is simulated for each frame separately as iid Poisson with intensity $\lambda _B=0.1$. For comparisons across different cases, the same simulated background frames are repeated in all examples.

• The dataset ${\mathbf{D }}$ is $64\times 64 \times 20$, a 3D grid of 20 frames of independent Poisson random variables totaling approximately 82K samples. For reference, the sample mean of the 20-frame background grid without hot spots is 0.0999 and the variance is 0.1001.

• To simulate a hot spot H, a subgrid is generated ${\mathbf{iid }}$ with Poisson intensity $\lambda _H$, $\lambda _H > \lambda _B$, and embedded in a background frame. Locations and shapes of hot spot subgrids within a $64\times 64$ frame and across the 20 frames are described below. To compare results when hot spots are added or deleted, the same simulated values for a hot spot are used regardless of whether other hot spots are present.

• The simulated hot spot geometries are a point source, a $2\times 2$ subgrid, and a larger angled shape with different Poisson intensities in two parts.

Hot spots $H_1$ and $H_2$

Two hot spots $H_1$ and $H_2$ are embedded in each frame of the 3D grid of background samples. $H_1$ is a point source—the same single grid point in each $64\times 64$ frame with iid sample values at Poisson intensity $\lambda _{H_1}=4.$ $H_2$ is a square $2\times 2$ subgrid with $\lambda _{H_2} = \lambda _{H_1}/4 = 1.$

Figure 1a is a 3D bar plot of the sample mean of the $64\times 64\times 20$ dataset ${\mathbf{D }}$ in its tensor mode-3, which is the $64\times 64$ matrix of sample mean values in the 20 frames. Figure 1b is the corresponding plot for the mode-3 mean values of tensor $\mathbf{M }$ that maximizes the Poisson likelihood $P({\mathbf{D }}|\mathbf{M })$.

Fig. 1

3D bar plots for two hot spots $H_1,H_2$. a Data tensor ${\mathbf{D }}$. b Computed tensor $\mathbf{M }.$

As noted in “Notation and computed tensor M” section, the sample means are equal and the sample variance in a mode of ${\mathbf{D }}$ is larger than in the same mode of $\mathbf{M }$. The sample mean is $\mu _{\mathbf{D }} = \mu _{\mathbf{M }}=0.1019$ and the sample variances are

$$\begin{aligned} \begin{array}{lc|ccc} &{}&{}{\mathbf{D }}&{}\mathbf{M }\\ \hline \text{ mode-3 }&{}&{}0.0103&{}0.0053\\ \text{ mode-2 }&{}&{}0.0018&{}0.0008\\ \text{ mode-1 }&{}&{}0.0017&{}0.0008\\ \end{array}\end{aligned}$$

Hot spots $H_1,H_2,H_3$

Figure 2 is the same layout as Fig. 1 but adds an additional hot spot $H_3$, a right-angled subgrid with its two arms having different Poisson intensities. The intensity in the arm parallel to the mode-1 axis is $\lambda _{H_1}/8 = 0.5$ and in the other arm is $\lambda _{H_1}/4= 1 =\lambda _{H_2}$.

Fig. 2

3D bar plots for three hot spots $H_1,H_2,H_3$. a Data tensor ${\mathbf{D }}$. b Computed tensor $\mathbf{M }.$

The sample mean is $\mu _{\mathbf{D }}=\mu _{\mathbf{M }}=0.1072$ and the sample variances are

$$\begin{aligned}\begin{array}{lc|ccc} &{}&{}{\mathbf{D }}&{}\mathbf{M }\\ \hline \text{ mode-3 }&{}&{}0.0150&{}0.0107\\ \text{ mode-2 }&{}&{}0.0021&{}0.0010\\ \text{ mode-1 }&{}&{}0.0022&{}0.0011\\ \end{array}\end{aligned}$$

Phase congruency in a 2D grid

The method of phase congruency

Both the data tensor ${\mathbf{D }}$ and the optimal Poisson-likelihood tensor $\mathbf{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, non-interactive) 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 ${\mathbf{D }}$ and $\mathbf{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\times Y$ matrix, the computation returns an $X\times 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.

Simulation

Two and three hot spots

Figure 3 shows phase congruency computed for the data ${\mathbf{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\times 64$ mode-3 sample mean values, the sample values averaged over the 20 frames. The 2D grid coordinates of the two closed contours in Fig. 3b define the location and outline of two potential hot spots.

Fig. 3

Phase congruency in data tensor ${\mathbf{D }}$ in Fig. 1. a, c Scatter plots for the 20 frames along mode-3. b 3-level contour plot

Let $N_H$ denote the sum of the sample values in a potential hot spot H. The integer $N_H$ is the Poisson count observed in the subgrid H of the larger grid ${\mathbf{D }}$. The counts enclosed by the respective contours of $H_1$ and $H_2$ in Fig. 3b are $N_{H_1} = 83$ (sample mean 4.15 in the 20 frames) and $N_{H_2} = 92$ (sample mean 1.15). The total count in ${\mathbf{D }}$ outside these two contours is 8010 (approximate background sample mean 0.098 in the 20 frames).

Phase congruency in the mode-1 and mode-2 projections is additional information. Figure 3a, c are scatter plots of phase congruency of the respective sample mean values (a) in the $20\times 64$ mode-2 projection retaining mode-3 and mode-1 axes and (c) in the $64\times 20$ mode-1 projection retaining mode-2 and mode-3 axes. These scatter plots are contour magnitudes with the extremes (lowest and highest) omitted. The mode-3 index steps from the newest frame 20 to the oldest frame 1, bottom-to-top in (a) and left-to-right in (c).

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 $\mathbf{M }$ in the same layout as Fig. 3. The two closed contours in Fig. 4b correlate strongly with Fig. 3b; however, the 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.

Fig. 4

Phase congruency in computed tensor $\mathbf{M }$ in Fig. 1. a, c Scatter plots for the 20 frames along mode-3. b 3-level contour plot

Figure 5 shows phase congruency in the data ${\mathbf{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 $\mathbf{M }$. $H_1$ is occluded by $H_3$ in the $20\times 64$ projection in (a), but the count in $H_1$ per frame accumulates in that mode-1 projection. The gap in scatter plot (a) at mode-1 indices 46 to 49 reflects the lower count in that arm of $H_3$ without $H_1$. The other arm of $H_3$ has Poisson intensity 1 instead of 0.5 and averages twice the count.

Fig. 5

Phase congruency in data tensor ${\mathbf{D }}$ in Fig. 2. a, c Scatter plots along mode-3. b 3-level contour plot

Fig. 6

Phase congruency in computed tensor $\mathbf{M }$ in Fig. 2. a, c Scatter plots along mode-3. b 3-level contour plot

The contours in Fig. 5b delineate the three subgrids $H_1, H_2, H_3$. Their respective counts are $N_{H_1} = 83$ (sample mean 4.15 in the 20 frames), $N_{H_2} = 92$ (sample mean 1.15), $N_{H_3}=489$ (sample mean 0.81). The total count in ${\mathbf{D }}$ outside these three subgrids is 7521 (approximate background sample mean 0.093).

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 ${\mathbf{D }}$ and $\mathbf{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 ${\mathbf{D }}$ in Fig. 7 has $H_2$ shielded in frames 1 to 9, then unshielded in frames 10 to 20. Figure 8 is the same layout for the computed tensor $\mathbf{M }.$ Qualitatively speaking, Figure 8 assists the localization of $H_2$ with better definition in the scatter plots than Fig. 7.

Fig. 7

Phase congruency in data tensor ${\mathbf{D }}$ with unshielding of $H_2$ at frame 10 (mode-3 index 10). a, c Scatter plots along mode-3. b 3-level contour plot

Fig. 8

Phase congruency in computed tensor $\mathbf{M }$ with unshielding of $H_2$ at frame 10 (mode-3 index 10). a, c Scatter plots along mode-3. b 3-level contour plot

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 ${\mathbf{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.

Fig. 9

Phase congruency in data tensor ${\mathbf{D }}$ with apparent motion of $H_2$ at frame 10 (mode-3 index 10). a, c Scatter plots along mode-3. b 3-level contour plot with dotted-line arrow in direction of motion

Fig. 10

Phase congruency in computed tensor $\mathbf{M }$ with apparent motion of $H_2$ at frame 10 (mode-3 index 10). a, c Scatter plots along mode-3. b 3-level contour plot with dotted-line arrow in direction of motion

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 ${\mathbf{D }}$ of 2D frames of Poisson count data, a 3D tensor factorization $\mathbf{M }$ is computed to maximize the conditional Poisson likelihood $P({\mathbf{D }}|\mathbf{M })$. The maximization of this likelihood is achieved by minimizing a Kullbach–Leibler divergence function [25, 40], specifically, $\mathbf{M }$ is computed in tensor math to minimize the function $\sum _i m_i - x_i \text{ log } m_i$ for index i over the nonnegative integers $x_i$ in ${\mathbf{D }}$ and the corresponding nonnegative reals $m_i$ in $\mathbf{M }$.

• 2. Phase congruency is computed for the two grids ${\mathbf{D }}$ and $\mathbf{M }$ projected in their three tensor modes. Phase congruency provides an invariant, normalized, numerical measure of “events” in ${\mathbf{D }}$ and in $\mathbf{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 $\mathbf{M }$ augments the data tensor ${\mathbf{D }}$ for analysis of phase congruency. There is strong correlation between ${\mathbf{D }}$ and $\mathbf{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 ${\mathbf{D }}$ due to Poisson clumping of background counts, even at a low Poisson rate. The computation of $\mathbf{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 ${\mathbf{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].

Appendix

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 ${\mathbf{D }}$ of independent Poisson samples.

In contemporary numerical math, an $n{\text{th}}$-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\times M_2$ matrix is an order-2 tensor. If tensor ${\mathbf{D }}$ is an $M_1\times M_2 \times \cdots \times 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 $n\text{{th}}$-order tensors [29, 33]. Further results in tensor math concern factorizations for least-squares error (LSE); specifically, given tensor ${\mathbf{D }}$, the tensor $\hat{\mathbf{D }}$ is computed in a designated factorized form to minimize the norm $||{\mathbf{D }} - \hat{\mathbf{D }}||_F$ [32, 33]. (The Frobenius norm $||\mathbf{Y }||_F$ of tensor $\mathbf{Y }$ is the square root of the sum $\sum _i y_i^2$ over all the values $y_i$ in $\mathbf{Y }$ [30]).

Implicit in using LSE for a rank-1 tensor factorization of ${\mathbf{D }}$ is an assumption of Gaussian variation in the data [24]. An alternative is to compute a Poisson model $\mathbf{M }$ for tensors of nonnegative counts which might be sparse [24, 25]. This includes nonnegative arrays of radiation data ${\mathbf{D }}$ which is Poisson in distribution.

Given order-n tensor ${\mathbf{D }}$ of nonnegative integer counts, the tensor $\mathbf{M }$ is computed to maximize the conditional Poisson likelihood $P({\mathbf{D }}|\mathbf{M })$. $\mathbf{M }$ is a tensor of nonnegative real numbers the same size as ${\mathbf{D }}$. Letting index i denote a 3D-location in ${\mathbf{D }}$ and $x_i$ [or $m_i$] the value in ${\mathbf{D }}$ [or $\mathbf{M }$] at i, then the likelihood $P({\mathbf{D }}|\mathbf{M })$ is maximized by minimizing the function

$$\sum _i m_i - x_i \text{ log } m_i,$$

the negative log of $P({\mathbf{D }}|\mathbf{M })$ except for an additive constant that depends only on the dataset ${\mathbf{D }}$. The set $\{m_i\}$ of local Poisson parameters in $\mathbf{M }$ optimizes $P({\mathbf{D }}|\mathbf{M })$ assuming independent samples of Poisson random variables in ${\mathbf{D }}$.

The iterative computation developed in [24, 25] computes the optimal tensor $\mathbf{M }$ in a rank-1 factorized form. For an $n\text{{th}}$-order tensor ${\mathbf{D }}$ and a user-specified number of components R, the method computes scalors $\alpha _r$ for $1 \le r \le R$ and matrices $A^{(i)}$ of respective sizes $M_i\times R$ for $1 \le i \le n$; then the tensor $\mathbf{M }$ is

$$\mathbf{M } = \sum _{r=1}^R \alpha _r a_r^{(1)}\circ a_r^{(2)}\circ \cdots \circ a_r^{(n)}$$

where $a_r^{(i)}$ is the $r\text{th}$ column of matrix $A^{(i)}$ and $\circ$ is outer product. The outer product of n matrices creates an $n\text{th}$-order tensor.

Convergence of the iteration to a tensor $\mathbf{M }$ maximizing $P({\mathbf{D }}|\mathbf{M })$ has been proved for mild conditions on the nonzero values in ${\mathbf{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 $\mathbf{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,22,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 ${\mathbf{D }}$ or from an array of Poisson parameters like computed tensor $\mathbf{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\times 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 $\lambda _{min}=3$ are suppressed.

The maximum and minimum moments of phase congruency covariance are $X\times 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 ${\mathbf{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.