# Identifying diffuse spatial structures in high-energy photon lists

Minjie Fan  Department of Statistics, University of California, Davis,
One Shields Avenue,
Davis, CA 95616, USA  Jue Wang  Department of Statistics, University of
California, Davis,
One Shields Avenue,
Davis, CA 95616, USA  Vinay L. Kashyap Center for Astrophysics $|$ Harvard &
Smithsonian,
60 Garden Street,
Cambridge, MA 02138, USA  Thomas C. M. Lee Department of Statistics,
University of California, Davis,
One Shields Avenue,
Davis, CA 95616, USA  David A. van Dyk Statistics Section, Department of
Mathematics, Imperial College London,
180 Queen’s Gate,
London, SW7 2AZ, UK  Andreas Zezas  Physics Department, deUniversity of Crete,
P. O. Box 2208, GF-710 03,
Heraklion, Crete, Greece

###### Abstract

Data from high-energy observations are usually obtained as lists of photon
events. A common analysis task for such data is to identify whether diffuse
emission exists, and to estimate its surface brightness, even in the presence
of point sources that may be superposed. We have developed a novel non-
parametric event list segmentation algorithm to divide up the field of view
into distinct emission components. We use photon location data directly,
without binning them into an image. We first construct a graph from the
Voronoi tessellation of the observed photon locations and then grow segments
using a new adaptation of seeded region growing, that we call Seeded Region
Growing on Graph, after which the overall method is named SRGonG. Starting
with a set of seed locations, this results in an over-segmented dataset, which
SRGonG then coalesces using a greedy algorithm where adjacent segments are
merged to minimize a model comparison statistic; we use the Bayesian
Information Criterion. Using SRGonG we are able to identify point-like and
diffuse extended sources in the data with equal facility. We validate SRGonG
using simulations, demonstrating that it is capable of discerning irregularly
shaped low surface-brightness emission structures as well as point-like
sources with strengths comparable to that seen in typical X-ray data. We
demonstrate SRGonG’s use on the Chandra data of the Antennae galaxies, and
show that it segments the complex structures appropriately.

Astrostatistics — Astrostatistics techniques — Astrostatistics strategies —
X-ray astronomy – Galaxy structure

††journal: AJ††facilities: Chandra††software: CIAO

##  1 Introduction

A challenge often encountered in high-energy astronomical analysis is that the
images are photon starved and sparse, and contain many ‘empty’ pixels. Unlike
photon-rich images encountered at longer wavelengths, complex features in
X-ray and $\gamma$-ray data are difficult to recognize, characterize, and
analyze. Working directly with Poisson distributed photon counts, while
simultaneously separating out the contribution of the background, is a
difficult process, especially when trying to detect faint non-uniform
emission, or separating faint point sources from larger scale diffuse
emission. Finding the boundaries of extended structures is thus a challenging
problem. Such complex structures are common in high-energy astronomical images
and include, for example, shock fronts, knots in supernova remnants, regions
of diffuse emission in galaxies, point sources embedded in diffuse emission or
conglomerates of point sources, entire galaxies or groups/clusters of
galaxies, jets, or star forming regions which appear to be extended in the
X-ray band even with intermediate resolution ($\lesssim 0.5\arcmin$) X-ray
telescopes.

The analysis of extended X-ray sources is critical for several areas of
astrophysics. The spatial scales of extended emission contain information
regarding the physical processes that lead to their formation, while their
boundaries are often determined by their physical environments. Therefore,
identifying the boundary of these regions in a data-driven, rather than a
model-driven, fashion is necessary for the scientifically valuable results. In
such cases, a primary goal of the researcher is to segment the image into
regions with similar properties and to analyze each segment individually.

Multiscale methods like wavelets (Starck et al., 2002) for point source
detection has been efficiently implemented for X-ray images (Freeman et al.,
2002), and matched-filter techniques have been successfully used to detect
galaxy clusters in ROSAT data (Vikhlinin et al., 1998), but extended
structures remain difficult to find and characterize in these low counts
Poisson data. Other techniques are generally optimized for high S/N images:
they apply adaptive binning, or set S/N thresholds to smoothed images with
point sources removed (e.g., Sanders & Fabian, 2001; Sanders, 2006); adapt
methods developed for the analysis of cosmic microwave background images
(e.g., Bobin et al., 2016); limit themselves to restrictive assumptions like
modeling a combination of point sources ((E)BASCS; Jones et al., 2015; Meyer
et al., 2021), or require spectral model similarity across the field of view
(Picquenot et al., 2019, 2021). Currently, most astronomical images with
complex structures that are processed for public display use some form of
flux-non-conserving adaptive smoothing (Ebeling et al., 2006). This approach
is inadequate for scientific analysis. Previous efforts at extended source
detection using Voronoi tessellation techniques have been limited by
computational cost and the imposition of global thresholding schemes (e.g.,
vtpdetect, Ebeling & Wiedenmann, 1993). Methods akin to seeded region growing
(cf. SrcExtractor, Bertin & Arnouts 1996; NoiseChisel, Akhlaghi & Ichikawa
2015) and Machine Learning (e.g., Morpheus, Hausen & Robertson 2020; Mask
R-CNN, Farias et al. 2020; galmask, Gondhalekar et al. 2022) have been used
for the identification of features in optical images of galaxies. However, the
Poisson nature and the sparsity of the X-ray data requires statistically
better targeted methods.

Here we develop a new method that combines aspects of Voronoi tessellation
with region growing by using neighbor similarity clustering. The method can be
applied to X-ray data and provides both separation between different
structures in a complex image and well-defined apertures to perform
photometry. We describe the statistical model that underlies the method in
Section 2, and specific implementation details including computational methods
in Section 3. We carry out several simulations to test the limits of
applicability of the algorithm in Section 4, and apply it to Chandra data of
the Antennae galaxies in Section 5. We discuss how and when the algorithm may
be best used in 6, and summarize our work in Section 7.

##  2 Statistical Methodology

We have developed a method that iteratively aggregates contiguous sets of
photons into distinct regions based on similarity of their surface brightness.
We employ a likelihood based method to obtain a piece-wise constant estimate
of the surface brightness across the image; the likelihood function is derived
in Section 2.1. The method starts with the high-resolution segmentation of the
spatial distribution of the events based on the Voronoi cells described in
Section 2.2 and combines segments by optimizing the Bayesian Information
Criterion (BIC) given in Section 2.3. Table 1 provides a glossary of our
notation.

Table 1: Glossary of variables and notation. Notation | Description
---|---
$\buildrel\rm iid\over{\sim}$ | Independent and identically distributed
$\hat{\zeta}$ | Estimate of a generic parameter $\zeta$
$\mathcal{F}$ | Field of view, a bounded domain in $\mathbb{R}^{2}$ that contains the observed photons
$n$ | Number of observed photons
$\mathbf{X}=\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\\}$ | Observed location of the $n$ photons, denoting their (sky) coordinates
$K$ | Number of segments within $\mathcal{F}$
$\mathbf{S}=\\{\mathcal{S}_{1},\ldots,{\mathcal{S}_{K}}\\}$ | Partition of $\mathcal{F}$, where $\mathcal{S}_{k}$ is the domain for segmented region $k$
${{\rm Num}(\mathcal{S}_{k})}$ | Number of photons observed in segment $\mathcal{S}_{k}$
${{\rm Area}(\mathcal{S}_{k})}$ | Area of segment $\mathcal{S}_{k}$
$m_{\rm seg}$ | Number of free parameters per segment
$\lambda(\mathbf{x})$ | Poisson intensity at location $\mathbf{x}$
$\bm{\mathbf{\lambda}}=\\{{\lambda_{1}},\ldots,{\lambda_{K}}\\}$ | Collection of the Poisson intensities over each of the segments
$\mathbf{\Lambda}$ | Total integrated intensity of $\lambda(\mathbf{x})$ over $\mathcal{F}$
$\mathcal{V}_{i}$ | Voronoi cell defined by photon $i$
${{\rm Area}(\mathcal{V}_{i})}$ | Area of Voronoi cell $\mathcal{V}_{i}$
$\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}$ | The Voronoi estimator of the Poisson intensity across Voronoi cell $i$, with $\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}=1/{{\rm Area}(\mathcal{V}_{i})}$
${\sl f}(\mathbf{X},n)$ | Joint probability mass/density function of the observed number of photons and their locations
$\mathcal{L}(K,\mathbf{S},\bm{\mathbf{\lambda}}\ |\ \mathbf{X})$ | Log-likelihood of the modela
$\mathcal{L}_{\rm profile}(K,\mathbf{S}\ |\ \mathbf{X})$ | Profile log-likelihooda obtained by replacing $\mathbf{\lambda}$ in $\mathcal{L}(K,\mathbf{S},\mathbf{\lambda}\ |\ \mathbf{x})$ with its estimate $\hat{\mathbf{\lambda}}$
${\rm BIC}(K)$ | Bayesian Information Criterion as a function of the number of segments, $K$
$m_{\rm grid}$ | In seed specification, the number of grid points used in a regularly spaced grid
$m_{\rm graph}$ | Number of photons included in initial subgraph for each seed
$m_{\rm nn}$ | Number of nearest neighbors over which local maxima are searched for to specify seeds
$m_{\rm R}$ | Number of strata used during Voronoi-area stratified sampling to specify seeds
$m_{\rm vorthr}$ | Minimum number of photons required for a Voronoi-area derived seed to be accepted
$a:$ Notation is reversed, by convention, for $\mathcal{L}(b|a)$ compared to
conditional probabilities; e.g.,
${\rm p}(a|b)$ represents the probability of $a$ given $b$.

###  2.1 Statistical Model

We consider an event list composed of $n$ photons observed in a bounded domain
that defines the field of view, $\mathcal{F}\subset\mathbb{R}^{2}$. Ignoring
instrumental pixelization, we model the set of sky coordinates for the $n$
photons,

$\mathbf{X}=\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\\}\,,$ (1)

via an inhomogeneous Poisson process with intensity function
$\lambda(\mathbf{x})\geq 0$. The intensity function must be integrable over
$\mathcal{F}$, i.e.,
$\mathbf{\Lambda}=\int_{\mathcal{F}}\lambda(\mathbf{x})d\mathbf{x}$ must be
finite. For simplicity, we assume that the intensity function is piece-wise
constant. Specifically, we assume we can partition $\mathcal{F}$ into $K$
segments, denoted $\mathbf{S}=(\mathcal{S}_{1},\ldots,\mathcal{S}_{K})$, where
the Poisson intensity is constant on each $\mathcal{S}_{k}$. Since
$\mathbf{S}$ partitions $\mathcal{F}$, the $\mathcal{S}_{k}$ together cover
$\mathcal{F}$ and each pair is disjoint. For a given set of non-negative
intensities, $\bm{\mathbf{\lambda}}=\\{{\lambda_{1}},\ldots,{\lambda_{K}}\\}$,
we can then express the intensity function as

$\lambda(\mathbf{x})=\sum_{k=1}^{K}{\lambda_{k}}\mathbbm{1}_{\mathcal{S}_{k}}(\mathbf{x})\,,$
(2)

where $\mathbbm{1}_{\mathcal{S}_{k}}(\mathbf{x})$ is an indicator function
that takes value 1 if $\mathbf{x}\in\mathcal{S}_{k}$ and is otherwise 0.

A property of the inhomogeneous Poisson process is that the number of photons,
${{\rm Num}(\mathcal{S}_{k})}$, recorded in segment $\mathcal{S}_{k}$ with
area ${{\rm Area}(\mathcal{S}_{k})}$ follows a Poisson distribution with mean

${{\rm
Area}(\mathcal{S}_{k})}\cdot{\lambda_{k}}=\int_{\mathcal{S}_{k}}\lambda(\mathbf{x})d\mathbf{x},~{}\hbox{
for }~{}k=1,\ldots,K\,,$ (3)

with ${\lambda_{k}}\geq 0$. Likewise, the total photon count is distributed
$n\sim{\rm Pois}(\mathbf{\Lambda})$. Another property is that, given $n$, the
sky coordinates, $\mathbf{x}_{i}$ are independent and identically distributed
(iid) with (normalized) probability density function
$\lambda(\mathbf{x})/\mathbf{\Lambda}$ (e.g., Chiu et al., 2013). This means
that the $\mathbf{x}_{i}$ are distinct – no two photons can have the same
recorded coordinates. (The discrete nature of detectors means that
occasionally, two photons are recorded with identical coordinates. In this
case, we add a very small random scatter, $\sim$10-6.) Our goal is to estimate
the number of segments, $K$, the segments, $\mathbf{S}$, and their respective
intensities, $\bm{\mathbf{\lambda}}$.

Thus far, we have not discussed sources or background. If the field of view
includes multiple separated sources, we expect the piece-wise constant
intensity function to capture the intensity peaks associated with point and
extended sources. Between these sources (or around a single source) is the
background region. If the background intensity is constant across
$\mathcal{F}$ and the source region(s) is/are isolated within the field of
view, we expect a single large segment representing the background to encircle
the source regions and to extend to the boundary of $\mathcal{F}$. If the
background intensity varies slowly, we might find several large segments that
together comprise the background region. In any case, there are segments
associated with background and with sources. We do not attempt to classify the
segments in this regard. Of course, if there is a single large low-intensity
segment encircling smaller higher intensity segments, it is easy enough to
identify the background with the large low-intensity segment.

We are particularly interested in the case where small-scale point-like
sources lie within a larger extended source as this is a challenging task for
existing methods. We do not distinguish between extended or point sources, and
in fact ignore the effect of telescope’s point spread function, assuming that
it is small compared to the size of the $\mathcal{F}$. Our method is agnostic
by design to the sizes of individual structures, and is thus capable of
isolating sources at all scales.

To derive the likelihood function, recall $n\sim{\rm Pois}(\mathbf{\Lambda})$
and given $n$ the $\mathbf{x}_{i}\buildrel\rm
iid\over{\sim}\lambda(\mathbf{x})/\mathbf{\Lambda}$. Thus, their joint
probability mass/density function under the inhomogeneous Poisson process is

$\displaystyle{\sl f}(\mathbf{x}_{1},\ldots,\mathbf{x}_{n},n)$
$\displaystyle=$ $\displaystyle{\sl f}(\mathbf{X}\mid n)\cdot{\sl f}(n)$ (4)
$\displaystyle=$
$\displaystyle\frac{1}{\mathbf{\Lambda}^{n}}\prod_{i=1}^{n}\lambda(\mathbf{x}_{i})\cdot\frac{\exp\left({-\mathbf{\Lambda}}\right)\mathbf{\Lambda}^{n}}{n!}$
$\displaystyle=$
$\displaystyle\frac{\exp\left({-\mathbf{\Lambda}}\right)}{n!}\prod_{i=1}^{n}\lambda(\mathbf{x}_{i})$

and their log-likelihood function is given by

$\mathcal{L}(K,\mathbf{S},\bm{\mathbf{\lambda}}\mid\mathbf{X},n)=\log{{\sl
f}(\mathbf{X},n)}\\\
=\sum\limits_{i=1}^{n}\log{\lambda(\mathbf{x}_{i})}-\int_{\mathcal{F}}\lambda(\mathbf{x})d\mathbf{x}-\log{n!}~{}~{},$
(5)

where we write out
$\mathbf{\Lambda}=\int_{\mathcal{F}}\lambda(\mathbf{x})d\mathbf{x}$. Replacing
$\lambda(\mathbf{x})$ by the piece-wise constant expression given in (2), we
have

$\mathcal{L}(K,\mathbf{S},\bm{\mathbf{\lambda}}\mid\mathbf{X},n)=\\\
\sum\limits_{\begin{subarray}{c}k=1\\\ {{\rm Num}(\mathcal{S}_{k})}\neq
0\end{subarray}}^{K}\kern-13.0pt{{\rm
Num}(\mathcal{S}_{k})}\log{{\lambda_{k}}}-\sum_{k=1}^{K}{{\rm
Area}(\mathcal{S}_{k})}{\lambda_{k}}-\log{n!}~{}~{}.$ (6)

Recall that ${{\rm Num}(\mathcal{S}_{k})}$ and ${{\rm Area}(\mathcal{S}_{k})}$
denote the number of photons in $\mathcal{S}_{k}$ and the area of
$\mathcal{S}_{k}$, respectively. (When ${{\rm Num}(\mathcal{S}_{k})}=0$, the
summand

${{\rm Num}(\mathcal{S}_{k})}\log{{\lambda_{k}}}$

is excluded from the first sum in (6).)

We aim to maximize $\mathcal{L}$ as a function of $K$, $\mathbf{S}$, and
$\bm{\mathbf{\lambda}}$ to obtain their maximum likelihood estimates. For
fixed $K$ and $\mathbf{S}$, $\mathcal{L}$ is maximized as a function of
$\bm{\mathbf{\lambda}}$ by

$\widehat{\lambda}_{k}={{\rm Num}(\mathcal{S}_{k})}/{{\rm
Area}(\mathcal{S}_{k})}~{}\hbox{ for }~{}k=1,\ldots,K\,.$ (7)

Plugging the $\widehat{\lambda}_{k}$ into $\mathcal{L}$, we obtain the profile
log-likelihood of $K$ and $\mathbf{S}$, i.e.,

$\mathcal{L}_{\rm profile}(K,\mathbf{S}\mid\mathbf{X},n)=\\\
\sum\limits_{\begin{subarray}{c}k=1\\\ {{\rm Num}(\mathcal{S}_{k})}\neq
0\end{subarray}}^{K}\kern-10.0pt{{\rm
Num}(\mathcal{S}_{k})}\log\left(\frac{{{\rm Num}(\mathcal{S}_{k})}}{{{\rm
Area}(\mathcal{S}_{k})}}\right)-n-\log{n!}.$ (8)

The same profile log-likelihood can be derived by modeling the data as a
mixture of uniform distributions. Allard & Fraley (1997) considered a special
case with a uniform background with a contiguous extended source superposed.

To estimate $\mathbf{S}$, we first deploy a (greedy) algorithm that finds an
optimal segmentation,
$\widehat{\mathbf{S}}(K)={\operatorname{arg\,max}}\mathcal{L}_{\rm
profile}(K,\mathbf{S})$, for fixed $K$, as described in Sections 2.2 and
refined in Section 3. In Section 2.3, we introduce a penalized version of
$\mathcal{L}_{\rm profile}$ that we maximize over $K$ to obtain final
estimates of the number of segments, $\widehat{K}$, and thereby of the
segments themselves, $\widehat{\mathbf{S}}(\widehat{K})$.

###  2.2 Estimating $\mathbf{S}$ via Voronoi Tessellation

Figure 1: Illustration of tessellation and triangulation of points. Left
Panel: 50 points drawn randomly from a multivariate normal
$N(\mathbf{0.5},\mathbf{0.2^{2}}{\cal{I}})$ are plotted as red dots, where
${\cal{I}}$ is a unit matrix. Middle Panel: The Voronoi cells generated by
these points are shown as cyan polygons surrounding each point. Right Panel:
The corresponding graph of Delaunay triangulation is shown as green line
segments connecting adjacent points.

Obtaining estimates of the segments, $\mathbf{S}$, requires us to constrain
the set of possible partitions. For any fixed $K$, for example, we can make
$\mathcal{L}_{\rm profile}$ arbitrarily large by including a segment that is
small enough to contain exactly one photon and shrinking the segment’s area
toward zero (since ${{\rm Area}(\mathcal{S}_{k})}$ appears in the denominator
of (8)). Similarly, any $\mathcal{S}_{k}$ with ${{\rm
Num}(\mathcal{S}_{k})}=0$ can have arbitrary shape since it does not
contribute to the profile likelihood. Since we cannot estimate the intensity
function at a higher resolution than the data, we only consider candidate
segments that include at least one photon.

One way to do constrain $\mathbf{S}$ is to only consider candidate segments
that consist of the Voronoi cells derived from the Voronoi tessellation of the
data, or the union of several Voronoi cells. The Voronoi tessellation of the
observed photons uniquely partitions $\mathcal{F}$ into $n$ convex cells,
denoted $\mathcal{V}_{i},~{}i=1,\ldots,n$, such that cell $\mathcal{V}_{i}$
contains exactly one photon, say $\mathbf{x}_{i}$, and consists of all
locations in $\mathcal{F}$ closer to photon $\mathbf{x}_{i}$ than to any other
photon. These cells are called Voronoi cells, and $\mathbf{x}_{i}$ is called
the nucleus of $\mathcal{V}_{i}$. Figure 1 gives an example of the Voronoi
tessellation of 50 photons drawn from a normal distribution truncated to the
unit square. The photon locations are plotted in the left panel and the
Voronoi cells in the middle panel. (We discuss the graph in the right panel in
Section 3.1.)

To avoid unclosed Voronoi cells near the border of the field of view, we
restrict the tessellation to Voronoi cells whose vertices are all in
$\mathcal{F}$. Based on the Voronoi tessellation, Barr & Schoenberg (2010)
introduced the Voronoi estimator
$\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}(\mathbf{x})=1/{{\rm
Area}(\mathcal{V}_{i})}$ for any location $\mathbf{x}\in\mathcal{V}_{i}$. They
show that under certain conditions, the Voronoi estimator is approximately
unbiased for the Poisson intensity $\lambda(\mathbf{x})$, and its sampling
distribution is approximately the inverse Gamma distribution111The probability
density function of an inverse Gamma distribution is
$\frac{b^{a}}{\Gamma(a)}x^{-a-1}\exp\left(-\frac{b}{x}\right)$, where $x>0$,
$a$ and $b$ are the shape and rate parameters, respectively, and
$\Gamma(\cdot)$ denotes the Gamma function. .

The algorithm that we propose to combine the Voronoi cells to form the
segments (by approximately maximizing $\mathcal{L}_{\rm profile}$ for each
fixed $K$) is detailed in Section 3. When $K$ is fixed in advance, this
algorithm can be used to estimate $\mathbf{S}$. To fit $K$ we use the method
in Section 2.3.

###  2.3 Estimating $K$ via the Bayesian Information Criterion

Unfortunately, the number of sources, $K$, cannot be reasonably estimated by
maximizing the profile likelihood, because (8) increases with $K$ and is thus
maximized by $K=n$, i.e., with the full set of Voronoi cells222Since $n$ is
fixed, maximizing $\mathcal{L}_{\rm profile}$ is equivalent to maximizing the
sum in (8), which can be written as
$\sum_{i=1}^{n}\log\hat{\lambda}(\mathbf{x}_{i})$ where
$\hat{\lambda}(\mathbf{x}_{i})$ is the local optimizer, ${{\rm
Num}(\mathcal{S})}/{{\rm Area}(\mathcal{S})}$, for the segment containing
$\mathbf{x}_{i}$. Increasing the number of segments allows for better local
optimization of local fluctuations and thus increases $\mathcal{L}_{\rm
profile}$. Of course, with too many segments, better fitting of local
fluctuations amounts to fitting noise, i.e., over-fitting the data.. We avoid
such over-fitting by adding a term to (8) that suitably penalizes model
complexity. Specifically, we use the so-called Bayesian Information Criterion
(BIC), which has been shown to produce statistically consistent results for
many model selection problems. For the current problem, the BIC is defined as

${\rm BIC}(K)=-2\mathcal{L}_{\rm
profile}(K,\widehat{\mathbf{S}}(K)\mid\mathbf{X},n)+Km_{\rm seg}\log n\,,$ (9)

where $m_{\rm seg}$ is the number of free/independent parameters per segment,
thus $m_{\rm seg}K$ is the total number of free parameters in the
model333Since the shape of the final segment is determined by the first $K-1$
segments, a more precise formulation of the total number of parameters in the
model is $m_{\rm seg}(K-1)+1$, where the intensity parameter of the final
segment is accounted for by the “+1”. The difference between this more precise
formulation and the one used in (9) is $(1-m_{\rm seg})n$, which does not
depend on any of the unknown parameters and thus does not affect estimation..
The BIC estimate for $K$ is given by

$\widehat{K}={\operatorname{arg\,min}}~{}{\rm BIC}(K)\,.$ (10)

For fixed $K$, optimizing BIC is equivalent to optimizing $\mathcal{L}_{\rm
profile}$, thus this estimate of $\mathbf{S}$ is equivalent to that described
in Section 2.2.

Unfortunately, because we are using a non-parametric model, $m_{\rm seg}$ is
not well-defined. Following Aue & Lee (2011), we set $m_{\rm seg}$ by
approximating the model by a parametric one using an assumed specific shape
for the segments. When the segments are close to ellipse-shaped, for example,
$m_{\rm seg}=6$ to account for the coordinates of the center, lengths of the
two axes, orientation, and intensity. When the segments are close to circular,
$m_{\rm seg}$ is reduced to $4$. Another possibility is to simply set $m_{\rm
seg}=1$ and the number of model parameters to the number of segments, which to
some extent reflects the overall model complexity (Magnussen et al., 2006),
but ignores the shapes of the segments. While these parametric approximations
allow us to assign a reasonable value to $m_{\rm seg}$, the model itself
remains non-parametric.

BIC is closely related to the “fitness function” used in the Bayesian block
method (Scargle et al., 2013), where model complexity is penalized via a
geometric prior on the number of sources, i.e., $p(K)=P_{0}\gamma^{K}$, with
$\gamma$ being a tuning parameter and $P_{0}$ a normalization constant.
Setting $\gamma=1/{n}^{Km_{\rm seg}}$ makes the fitness function equivalent to
the penalty term in BIC.

##  3 Algorithm for Combining Voronoi Cells into Segments

###  3.1 SRGonG: Seeded Region Growing on Graph

Using the Voronoi cells as building blocks, we start by proposing an algorithm
to estimate the segments in $\mathbf{S}$ for fixed $K$. The first step is to
identify pairs or groups of Voronoi cells that can potentially be combined. We
accomplish this via the dual graph of the Voronoi tessellation, known as the
Delaunay triangulation. This graph’s vertices are the centers of the Voronoi
cells, i.e., the photons, and its edges connect pairs of the adjacent Voronoi
cells. The right panel of Figure 1 depicts the graph derived from the Voronoi
cells in the middle panel. We assign vertex $\mathbf{x}_{i}$ the value of the
Voronoi estimator, denoted by
$\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}$, i.e., an estimate of
the intensity in Voronoi cell $\mathcal{V}_{i}$. Using the graph constructed
by the Delaunay triangulation, the problem of estimating the $\mathcal{S}_{k}$
can be naturally translated to the problem of graph segmentation, i.e.,
partitioning the graph into subgraphs such that the Voronoi cells therein form
a single segment, $\mathcal{S}_{k}$. Thus, each subgraph/$\mathcal{S}_{k}$ is
formed by a set of vertices/photons connected by a collection of edges in the
full graph. Since we assume that the intensity function is piece-wise
constant, all the vertices in each subgraph should share similar values of
$\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}$.

Unfortunately, finding $\mathbf{S}$ to maximize (8) for fixed $K$ remains an
intractable combinatorial optimization problem even when confined to
combinations of Voronoi cells. A distinct advantage of representing the
problem as graph segmentation is that we implicitly impose an additional
constraint that each $\mathcal{S}_{k}$ is a subgraph. In this way, traditional
image segmentation444Image segmentation is the process of separating an image
into a number of regions such that each region is composed of connected pixels
with similar characteristics, such as similar pixel values. methods can be
adapted and used to segment the graph. In particular, we propose the Seeded
Region Growing on Graph (SRGonG) method, which is similar to the original
Seeded Region Growing (SRG) method used for images except that the concept of
“neighbors” is determined by the edges of the graph instead of neighboring
pixels. The original SRG is proposed in Adams & Bischof (1994) and is extended
to several variants to deal with more complicated cases in Fan & Lee (2014).
SRGonG starts by identifying, either manually or automatically, a set of
initial seeds from the graph. Each seed can be a single vertex/photon or a
seeding subgraph, i.e., a set of connected vertices/photons.

For the moment, we present a simplified version of SRGonG that requires a
perfect set of seeds, i.e., a set with exactly one seed in each
$\hat{\mathcal{S}}_{k}$. Recall that $K$ is fixed, thus initially we assume
$K$ seeds. The details of seed specification in more realistic settings are
described in Section 3.2 and the full version of SRGonG (which requires an
extra step to merge segments and estimate $K$) is detailed in Section 3.3.

SRGonG grows the seeds into subgraphs by successively adding neighboring
vertices to them. More specifically, at each iteration, the method selects a
pair that consists of a growing subgraph, $S$, and one of its unassigned
neighboring vertices, $i$, such that

$\delta(i,S)=\left\lvert\log\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}-\log\\{{{\rm
Num}(\mathcal{S})}/{{\rm Area}(\mathcal{S})}\\}\right\rvert$ (11)

is minimized. This criterion compares the logarithm of the estimated
intensities of the subgraph and the neighboring vertex because
$\mathcal{L}_{\rm profile}$, which we aim to optimize, combines the segment-
specific intensity estimates on the log scale. The vertex in the pair with the
smallest difference is added to the corresponding subgraph. This process
finishes when all the vertices of the full graph are assigned to exactly one
subgraph. The Voronoi cells contained in the subgraphs give the final
segmentation of $\mathcal{F}$, i.e., $\widehat{\mathbf{S}}(K)$, for pre-
specified $K$.

In practice, we save the index, $i$, of the neighboring Voronoi cell that
minimizes $\delta(i,S)$ for each growing subgraph, $\mathcal{S}_{k}$, at each
iteration. This reduces the time complexity of the method to be linear in
terms of the number of photons.

###  3.2 Seed Specification

Since SRGonG begins by building out regions starting from a specified set of
seeds, the number and location of the seeds are important considerations. As
discussed in Section 3.1, we would ideally have exactly one seed within each
$\hat{\mathcal{S}}_{k}$. Unfortunately, this is not feasible in practice. A
brute force solution is to over-specify the seed set to the extreme, by
setting every photon location to be a seed, and devising an algorithm to merge
the seeds into segments. But merging such a large seed set would be
challenging in terms of both computational speed and statistical accuracy (see
discussion in Section 6.1). If the field being analyzed is known to have a
large number of point sources, or if the scientific question requires focusing
on point sources, then running a source detection algorithm first to find all
such sources and specifying all of them as seeds will be helpful. Here we
describe three generic strategies to specify smaller initial seed sets. These
strategies still over-specify the set in that they use a larger number of
seeds than the expected number of segments (but less so than setting each
photon to be a seed). Thus, after growing the seeds into subgraphs as
described in Section 3.1, we require a method to merge the resulting subgraphs
into segments; we describe our merging algorithm in Section 3.3.

#### Regular grid:

This method starts by overlaying a regular grid of $m_{\rm grid}$ points onto
the field of view, $\mathcal{F}$. For each grid point, we specify a seeding
subgraph composed of the $m_{\rm graph}$ photons closest to the grid point (in
terms of the Euclidean distance). We typically set $m_{\rm grid}$ and $m_{\rm
graph}$ so that their product is much smaller than $n$ to enable the seeded
regions to grow. Conflicts in the allocation of photons to seeding subgraphs
(e.g., when a single photon is among the $m_{\rm graph}$ closest to two or
more grid points) are broken by the order of assignment. The number of the
photons, $m_{\rm graph}$, assigned to each seeding subgraph can be increased
to stabilize the initial estimates of the growing subgraph intensities,
especially when the contrast (i.e., the ratio between the intensities of an
extended source and the background) is low. In practice, there is no
universally best choice for $m_{\rm grid}$ and $m_{\rm graph}$, as their
optimal values depend on factors including the number of observed photons $n$
and the complexity and number of true astronomical sources. To ensure a
sufficient number of photons for the seeding subgraphs, we require

$m_{\rm graph}{\leq}\frac{n}{m_{\rm grid}}\,.$

Since it is possible for a seed to fall on the boundary between two
distinguishable segments and adversely affect subsequent processing, we
propose an additional seed-rejection step. Specifically, we compare the range
of the Voronoi areas for each photon ${{\rm
Area}(\mathcal{V}_{k})}|_{\\{k=1,\ldots,{m_{\rm graph}}\\}}$ of a seeding
subgraph of size $m_{\rm graph}$ with the expected empirical $2\sigma$
confidence interval for a homogeneous distribution of photons (Møller, 1994,
Chapter 4.2), i.e.,

$\frac{1}{\widehat{\mathbf{\lambda}}_{s}}\pm
2\times\frac{0.53}{\widehat{\mathbf{\lambda}}_{s}}$ (12)

(Møller, 1994, Chapter 4.2), where

$\frac{1}{\widehat{\mathbf{\lambda}}_{s}}=\frac{1}{m_{\rm
graph}}\sum_{k=1}^{m_{\rm graph}}{{\rm Area}(\mathcal{V}_{k})}$ (13)

is the average of Voronoi areas for the photons in the seeding subgraph. Thus,
if the actual range of Voronoi areas ${{\rm Area}(\mathcal{V}_{k})}$ exceeds
the expected empirical confidence interval, the seeding subgraph is rejected.

#### Grid supplemented by local maxima:

If the regular grid used to generate the seeds is too sparse, some image
structures may not be captured in the segmentation. For example, if there is
not a grid point sufficiently near a point or extended source, the source may
be merged into the background or another source. One remedy is to include
additional seeds near the likely locations of sources. Sources induce an
elevated intensity over small spatial scales. Thus, locations of high photon
density are likely associated with sources. We propose to identify vertices
that are local maxima of the graph constructed by the Delaunay triangulation,
in the sense that the vertex value (i.e.,
$\widehat{\lambda}_{i}^{{\rm\mathcal{V}}_{i}}$) is greater than or equal to
that of its closest $k$ vertices (including itself), where closeness is
measured by the Euclidean distance. For each local maxima we find in this way,
we include a seeding subgraph composed of its closest $m_{\rm graph}$ vertices
(including itself).

#### Voronoi-area stratified sampling:

More complex schemes, designed to locate seeds over a broader range of surface
brightness, can also be devised. Methods such as Otsu’s thresholding (Otsu,
1979) can also be used to specify seeds or seeding subgraphs for point-like or
localized extended sources. As we discuss in Section 6.1, the grid
supplemented by local maxima is adequate to identify structures that exist at
a large variety of scales in astronomical data. Here, as an example case, we
describe a third method, which selects seeds via stratified sampling of the
photons, with strata determined by the areas of the Voronoi cells, ${{\rm
Area}(\mathcal{V}_{i})}$. Specifically, the photons are divided into $m_{\rm
R}$ strata bounded by equally spaced quantiles of the distribution of ${{\rm
Area}(\mathcal{V}_{i})}$. The number of strata depends on the sample size, but
we typically use $m_{\rm R}\approx{10-20}$. Clumps of near-neighbor photons
within each stratum are put together (see Appendix A) into a set of labeled
groups such that spatially nearby photons within a given stratum are all
assigned the same label. If a given label is assigned to fewer than $m_{\rm
vorthr}$ photons (typically $m_{\rm vorthr}=5$), then the photons with this
labels are discarded for the purpose of seed specification; otherwise the
central photon555Consider the set, $L$, of photons with a given label. For
each photon $l\in L$, we calculate $d_{l}=\sqrt{{{\rm
Area}(\mathcal{V}_{l})}}+\sum_{k{\in}L}~{}d_{lk}$ where $d_{lk}$ is the
Euclidean distance between photons $l$ and $k$. That photon in $L$ with the
smallest $\\{d_{l}\\}$ is flagged as the central photon among the photons in
$L$. This measure of centrality is better than computing a centroid as it
ensures that the seed is guaranteed to be included inside the labeled region
even when the region shape is complex, and that the seed is unambiguously
assigned to one of the photons. amongst those with each label is set as a
seed. Subgraphs are then constructed for each retained seed photon in the same
manner as described in Section 3.1.

Data: Coordinates of observed photons
$\mathbf{x}_{i}=(x_{1i},x_{2i}),i=1,\ldots,n$ in field of view $\mathcal{F}$.

Result: Piece-wise constant estimate of intensity function with a segmentation
of $\mathcal{F}$ into regions of constant intensity,
$\widehat{\mathbf{S}}=(\hat{\mathcal{S}}_{0},\ldots,\hat{\mathcal{S}}_{\widehat{K}})$.

begin

[0pt]Use Voronoi tessellation to obtain a graph whose vertices are the
observed photons with the Voronoi estimators
$\widehat{\lambda}_{i}^{{\rm\mathcal{V}}_{i}}$ as their values.

[0pt]Using a method in Section 3.2, specify the initial seeds for subgraph
growing.

[0pt]Grow seeds into subgraphs that over-segment the entire graph:

while  _there are unassigned vertices_ do

Select a pair of a growing subgraph $S$ and one of its neighboring vertices
$i$ such that
$\delta(i,S)=\left\lvert\log\widehat{\mathbf{\lambda}}_{i}^{{\rm\mathcal{V}}_{i}}-\log\\{{{\rm
Num}(\mathcal{S})}/{{\rm Area}(\mathcal{S})}\\}\right\rvert$ is minimized.

Add the vertex in the pair with the smallest difference to the corresponding
subgraph.

[0pt]Greedily merge subgraphs by minimizing BIC at each merger to obtain a
nested sequence of segmentations.

[0pt]Finally, set $\widehat{K}$ and $\widehat{\mathbf{S}}$ to the values of
the nesting level with the smallest BIC. $\widehat{K}$ is the final
segmentation of $\mathcal{F}$.

Algorithm 1 Seeded Region Growing on Graph (SRGonG)

###  3.3 Subgraph Merging

Using one of the seed sets of Section 3.2 to grow subgraphs as described in
Section 3.1 leads to an over-segmented graph since the number of seeds is
invariably more than the predetermined $K$ or the $\widehat{K}$ that optimizes
BIC. To merge the subgraphs into segments, we propose a subgraph merging
method that aims to minimize BIC. Similar ideas were used by Lee (2000) and
Peng et al. (2011) in image segmentation. Specifically, the subgraph merging
method starts by computing the BIC for the over-segmented graph and then
iteratively selects two neighboring subgraphs, merges them, and recomputes
BIC. The two merged subgraphs are selected so that their merger gives the
largest decrease (or the smallest increase) of BIC among all possible merges
(of neighboring subgraphs). In this sense, this is a greedy algorithm. We
continue merging the subgraph until all subgraphs are merged into the entire
graph, except when $K$ is fixed, in which case we stop when $K$ segments
remain. In this way, we obtain a sequence of nested segmentations,
$\\{\widehat{\mathbf{S}}(K),K=1,\ldots,n\\}$, each with a BIC value. Finally,
we set $\widehat{K}$ and $\widehat{\mathbf{S}}$ to the values of the nested
level with the smallest BIC.

At each iteration, we use an updating formula to speed the computation of BIC.
Consider the graph segmentation associated with $K$ and the graph segmentation
after merging two of the subgraphs of $\mathbf{S}(K)$ and label the merged
subgraphs $i$ and $j$, with $1\leq i<j\leq K$. This merger decreases BIC by

$\displaystyle\Delta\mbox{BIC}_{K,i,j}$ $\displaystyle=$
$\displaystyle\mbox{BIC}(K)-\mbox{BIC}({K-1})$ (14) $\displaystyle=$
$\displaystyle 2\ {{\rm Num}(\mathcal{S}_{i})}\log\frac{{{\rm
Num}(\mathcal{S}_{i{\cup}j})}{{\rm Area}(\mathcal{S}_{i})}}{{{\rm
Area}(\mathcal{S}_{i{\cup}j})}{{\rm Num}(\mathcal{S}_{i})}}$ $\displaystyle+2\
{{\rm Num}(\mathcal{S}_{j})}\log\frac{{{\rm Num}(\mathcal{S}_{i{\cup}j})}{{\rm
Area}(\mathcal{S}_{j})}}{{{\rm Area}(\mathcal{S}_{i{\cup}j})}{{\rm
Num}(\mathcal{S}_{j})}}$ $\displaystyle+m_{\rm seg}\log n\,,$

where $i{\cup}j$ denotes the union of photons that belong to subgraphs $i$ and
$j$. The complete procedure for SRGonG is summarized in Algorithm 1.

Figure 2: Simulation studies with various shapes of the extended source (with
$\beta=1$ and $\sigma=30$; see Section 4). Each row corresponds to one of the
three image shapes (circular at top, polygonal zig-zag at middle, and semi-
annular at bottom) considered. Left Column: Simulated photon locations for one
simulation instance. Middle Column: The graph for the exemplar (first)
simulation of the left column is shown (grey lines), and the true segment
boundary is overlaid (solid green lines). The initial seeds are marked for the
$5{\times}5$ regular grid (black squares), local maxima (blue diamonds), and
rejected seeds (red triangles). Right Column: Segment boundaries from 10
additional simulations (grey solid lines) are overlaid on the photons from the
first simulation, which are marked as red for the extended source, and blue
for the point-like sources.  Figure 3: Demonstration of recovery of simulation
parameters in segmentation. The brightness of each segment recovered by SRGonG
for all of the simulations shown in Figure 2 is shown as the black histogram.
The expected brightness for each of the components (background, extended,
point-like) and shapes of the extended source (circle, polygonal zigzag, and
semi-annular) are marked as vertical lines with different line styles as
labeled. The background and extended emission brightnesses are recovered well,
but because of the large contrast, the brightness of point-like objects shows
a bias (see text).

m

m

Figure 4: Simulation studies for the various shapes of the extended source,
different exposure times, $\beta$, and different contrast ratios, $\sigma$.
Each point source has $\beta\cdot\sigma$ counts in a circle of radius 0.025,
each diffuse extended source has $10\cdot\beta\cdot\sigma$ counts spread
uniformly over an area $\approx$0.2, and there are $1000\cdot\beta$ counts in
the background, spread uniformly across $\mathcal{F}$. The three rows
correspond to the three extended source shapes in the rows of Figure 2. Left
Column: Boxplots of the adjusted Rand indices for $500$ replicates under each
simulation setting; Right Column: Histograms of the fitted number of sources
$\widehat{K}$ for the 500 replicates under each simulation setting, where the
true number of sources is highlighted in red.

##  4 Simulation studies

Our simulation study is conducted assuming a hypothetical instrument that
produces fields of view, $\mathcal{F}$, with two-dimensional coordinates on
the unit square. The simulations are designed to assess the performance of
SRGonG when applied to fields of view of point-like sources embedded in
extended sources of different shapes, while varying the exposure time (or
equivalently, the overall counts in the field) and the contrast between the
different components. In all our simulation settings, “point-like sources” are
circular sources of radius 0.025 and extended sources are of area $\approx$0.2
relative to the $\mathcal{F}$. We consider three “true images”: (a) four
point-like sources embedded within a circular extended source of radius 0.25
(covering an area $0.196$ of the unit square), (b) three point-like sources
embedded within a polygonal zig-zag shaped extended source comprising five
squares of size $0.2{\times}0.2$ (total area of $0.2$), and (c) three point-
like sources embedded within an arc-shaped extended source (a half-annular
shape with inner radius 0.2 and outer radius 0.4, and total area $0.189$). We
consider these three settings because point-like sources embedded within a
complex extended source are commonly observed in astrophysical fields of view,
as illustrated in Section 5 and the extended sources mimic typical
astronomical shapes. Furthermore, the variety of shapes and the contrasts
considered are a stringent test of the algorithm. Letting $\beta$ denote the
exposure time (in arbitrary units) and $\sigma$ denote the contrast level
between the different components, for each simulated $\mathcal{F}$ we
generated $\beta\sigma$ counts for each point-like source, $10\beta\sigma$
counts for the extended source, and $1000\beta$ counts in expectation for the
background, with the photons corresponding to each component distributed
uniformly over the area allocated to it. In Figure 2, we have adopted
$\beta=1$ and $\sigma=30$, so in all cases, the point-like sources have 30
photons, the extended shapes have 300 photons, and the background has
$\sim$Poisson(1000) photons, all distributed uniformly over their allocated
areas, with $\approx$200 background counts under the area of the extended
source. The contrast in surface brightness between the extended source and the
background is thus $\approx$1.5$\times$, which is sufficiently large on the
scale of the extended sources that the presence of the extended sources are
clearly recognizable. However, it is clear from inspection of Figure 2 that
local fluctuations can be sufficiently large as to make estimating the
boundary of the extended sources challenging.

The photons randomly generated for each of the settings are shown for one case
in the left column of Figure 2. The middle column shows shapes of the extended
sources are also shown overlaid on the corresponding Delauney triangulation
for each of the photons, as well as the seeds chosen for that case. The right
column shows the segmentation, with points colored blue and extended source
colored red, for the simulation in the left column, and superposed in grey
lines, the result of segmentations from 10 additional simulations. The
superpositions of the segment boundary lines over the expected lines of the
shapes of both the point-like and the extended sources show that while
fluctuations are present in individual simulations, on average the boundaries
are picked out well. A detailed examination of the locations of the boundaries
and their uncertainties requires modeling the boundaries, and we defer
discussion to follow-up work (J. Wang, et al., in preparation). Here, we
demonstrate that the components are well recovered in all cases. We show the
distribution of the segment brightnesses found for all the simulations for all
three cases in Figure 3: the components are clearly separated, with
uncertainties of $\approx$10-15% on the expected brightness in each component.
We find that the brightness of the point-like component suffers from a bias
because of the tendency of the segment areas to preferentially encroach on the
much larger area of the surrounding extended source, thus causing a downward
shift in the estimated brightness.

In our simulation design, in addition to varying the three “true images”, we
also vary the exposure time with $\beta$ taking values 0.5, 1 and 2, and the
contrast with $\sigma$ taking values 10, 20 and 30. We simulate 500 fields of
view under each of the 27 resulting simulation settings666The number of source
counts was held fixed in all simulations, while the number of background
counts was generated as a Poisson with mean $1000\beta$ in order to explore
the effect of background fluctuations. Thus, the total counts in a given
dataset is $\beta\sigma+10\beta\sigma+{\rm Poisson}(1000\beta)$.. Each of the
13,500 simulated fields of view is analyzed with SRGonG, with initial seeds
specified following the “grid supplemented by local maxima” method of Section
3.2. The regular grid used for seed specification is 5-by-5, with seed size
$m_{\rm graph}=5$ and a neighborhood size of $k=50$ for finding local maxima.
Since the sources we are considering are simple, we set the BIC parameter to
be $m_{\rm seg}=4$; when more complicated shapes are expected, larger values
of $m_{\rm seg}$ should be used.

The second and third columns of Figure 2 show the initial seed specifications
and segmentation results for the first of the 500 fields of view generated
with $\beta=1$ and $\sigma=30$. All the point-like sources are clearly
identified. The fitted boundaries of the extended sources are generally quite
good, except for some mild leakage for the arc-shaped extended source. In
general, we expect sources with longer perimeters per unit area777A standard
measure of shape irregularity is the “perimeter index” of Angel et al. (2010)
which is defined to be the perimeter of a circle of area equal to that of the
shape divided by the actual perimeter of the shape. to be more challenging.
This is because photons nearer the boundary of a segment are more likely to be
misclassified than are those nearer the middle. Thus, more irregularly shaped
sources, such as the arc-shaped source in this simulation, are more
challenging, including for SRGonG.

Furthermore, several of the initial seeds placed by the regular grid happen to
fall near the boundary of arc-shaped source, which can also jeopardize the
performance of seed-based methods.

We use a clustering verification metric, specifically the _adjusted Rand
index_ (ARI; Hubert & Arabie, 1985), to assess the quality of the SRGonG
segmentations. The Rand Index (RI; Rand, 1971) quantifies how well a given
segmentation matches the ground truth segmentation. Specifically, each pair of
photons is classified as either (a) being in the same fitted segment and in
the same ground truth segment, (b) being in different fitted segments and in
different ground truth segments, or (c) not being in class (a) or (b). (For
the ground truth, the segments are the background, extended source, and each
point-like source.) The Rand Index is defined to be the number of photons
pairs in class (a) or (b), relative to the total number of photon pairs. Thus,
a perfect match to the ground truth results in RI=$1$. The Adjusted RI
corrects the RI such that accidental overlaps of segments due to chance are
accounted for, yielding values in the range $-1<\hbox{ARI}<+1$.

Figure 4 summarizes the ARI and the fitted value for the number of segments,
$\widehat{K}$ for the $500$ replicates under each of the 27 simulation
setting. For each of the three true images, as expected, the SRGonG
segmentation improves as either the exposure time or the contrast between the
brightness of the components increases. This is seen in the progression from
the top left panel to the bottom right panel of the right column of plots in
Figure 4: the method fails to identify the embedded point sources when there
are $\approx$5 counts in each source, but correctly identifies all components
in $\geq$70% of the cases (80% for the circular and polygonal cases) when
there are 60 counts in each point source. Similarly, the ARI increases to
close to one (i.e., perfect agreement between ground truth and SRGonG
segmentation, where the fitted number of segments equals the true number of
sources) as $\beta$ and $\sigma$ increase.

##  5 Application to Antennae Galaxies

Figure 5: Representations of Antennae data. Top Left: Scatterplot of the
photons in from six Chandra observations of the Antennae galaxies carried out
between December 2001 and November 2002. The x-axis is normalized to range
from 0 to 1, and the y-axis is normalized accordingly, by the same ratio, such
that the data are depicted in the form that is input to SRGonG. Top Right: The
Delaunay triangulation of the photons is shown, with each vertex marked by a
point which is color coded by the brightness determined as the reciprocal of
the Voronoi area (see scale at right). Bottom Left: The seeding subgraphs for
SRGonG via a regular grid, marked with black squares, supplemented by local
maxima marked with blue diamonds. Seeds discarded due to a strong indication
of being on the boundary of two segments are marked as red triangles. Each
subgraph is represented by one photon in the middle of the seeding subgraph.
Bottom Right: The SRGonG segmentation of the Antennae galaxies data, with the
red curves depicting the boundaries of the fitted segments. Three smaller
fields of view are highlighted by black boxes, labeled Region 1-3, and
magnified to show details in Figure 7.

The Chandra observations of the Antennae galaxies provide a good test case for
application of SRGonG. The X-ray data (see Figure 5, top left panel) shows
complex structures. Specifically, the data reveal several point sources and
extended regions, with several clumps of diffuse emission of different extent
and surface brightness, along with a population of unresolved point-like
sources superposed. Some of the point sources lie within the extended sources
(e.g., in the extended region at the bottom of the image) and some of the
extended sources are entangled with each other.

As a conservative scenario we used the first Chandra observation of the
Antennae galaxies obtained on December 1st 1999 (OBSID 315; Fabbiano et al.,
2001). The observation was performed with the ACIS-S detector for a total
exposure of 72 ksec. We process and screen the data (e.g. initial
calibrations, removal of strong background flares, selection of good grades)
as in Zezas et al. (2022; CIAO v3.2, CALDB v2.11). Again as a conservative
scenario, we use the full dataset without any screening for events of very low
or high energies which are dominated by background. The final dataset we use
consists of $\approx$50,700 events within a
$\sim{3.45}\arcmin\times{3.45}\arcmin$ region around the galaxy (screening for
events in the generally used 0.5-8.0 keV band would result in a reduction of
$\sim 43\%$ in the total number of counts).

Figures 5 show different depictions of these data, with the coordinates scaled
linearly to the range $[0,1]$, as is assumed by our implementation of SRGonG,
and processed to show the resulting Voronoi tessellation.

We apply SRGonG to these data in order to obtain statistically meaningful non-
parametric segmentations of the different clumps of diffuse emission, as well
as to separate diffuse and point-like emission sources. We apply the Voronoi
tessellation to the photons and construct the graph of Delaunay triangulation
(see top right panel of Figure 5). We specify the initial seeds for SRGonG via
a regular grid supplemented by local maxima (see Section 3.2; shown in bottom
left panel of Figure 5). We start with a regular 9$\times$9 grid (i.e.,
$m_{\rm grid}=81$), the initial estimates of which are stabilized by assigning
the $m_{\rm graph}=20$ nearest photons to each seed; these cover the large
scale variations in the data. Local maxima are determined over a neighborhood
size of $k=100$. The 419 seeds that result from this process is a sufficiently
large number to ensure that there is generally at least one seed in each
point-like or extended source or the background. Since we expect the segments
of the extended sources to be more irregularly shaped in the real data than in
the simulation, we choose a larger value of the BIC parameter, $m_{\rm seg}=6$
(see Equation 9); this corresponds to assuming that each segment has the
complexity of an ellipse.

The results of SRGonG are shown in the bottom right panel of Figure 5 (the
regions outlined in black are discussed in more detail in Section 6.1),
showing the boundaries of the fitted segments as thin red lines around the
black dots depicting the photons. SRGonG correctly segments areas with similar
surface brightness such that photons that correspond to these diffuse
components are grouped together. The photons that belong to each of these
segments can be trivially collected together for further analysis, depending
on the scientific question being explored. For instance, in Figures 6, we show
segment-wise maps of the fractional hardness ratios HR${}_{SM}=(S-M)/(S+M)$
and HR${}_{MH}=(M-H)/(M+H)$, where $S,M,H$ are counts in PI channels [35:61]
(${\approx}$0.5:0.9 keV), [62:82] (${\approx}$0.9:1.2 keV), and [83:135]
(${\approx}$1.2:2 keV), respectively. Notice that the maps clearly demonstrate
that the diffuse emission in the Antennae generally have softer spectra than
the point sources. Maps such as these can be used to identify the extent of
dust lanes in the Antennae system; e.g., the segments at $\sim$(0.4,0.22),
$\sim$(0.3,0.25), and $\sim$(0.6,0.75), which are characterized by harder
spectra than the surrounding segments, a characteristic of increased
absorption (cf., Zezas et al., 2006). Furthermore, notice that the southern
region (around Region 3 in the bottom left panel of Figure 5) is surrounded by
a halo of relatively soft X-ray emission, in agreement with the spectral
analysis of Baldi et al. (2006) who find emission from soft $\sim$0.6 keV
thermal emitting gas.

Figure 6: Fractional hardness ratios of counts in each of the SRGonG-
determined segments of the Antennae data. The images show
HR${}_{SM}=\frac{S-M}{S+M}$ (left) and HR${}_{MH}=\frac{M-H}{M+H}$ (right),
where $S,M,H$ are counts in the passbands $0.5-0.9$ keV, $0.9-1.2$ keV, and
$1.2-2$ keV, respectively. In both figures, redder colors indicate softer
spectra. The background is rendered in reddish pink in HRSM and light green in
HRMH.

##  6 Discussion

###  6.1 Performance on the Antennae data

Figure 7: Magnified views of the three regions of the Antennae galaxies
highlighted in the bottom right panel of Figure 5. The grey dots represent
individual photons, and the thin red curves represent the SRGonG segmentation
(see Section 5). Some of the larger segments are labeled with blue letters
(see text).

Here we discuss the quality of the SRGonG segmentation of the Antennae in
greater detail. To begin with, we note that SRGonG successfully identifies a
number of point-like sources, characterized by the presence of a large number
of photons within a small space. Several of these point-like sources are
superposed on extended diffuse emission and surrounded by complex structures.
Furthermore, unlike the case usually with methods that use piece-wise constant
models, the point-like sources are invariably defined by single segments and
not several concentric rings that approximate the typical profile of the PSF
where intensity increases from the wings inward to rise to a peak at the core.
Such cases are not entirely absent, however; see Figure 7, specifically the
sources at $\sim$(0.23,0.57), $\sim$(0.43,0.54) and $\sim$(0.31,0.46) in
Region 1, $\sim$(0.67,0.4) in Region 2, and $\sim$(0.49,0.12) and
$\sim$(0.37,0.07) in Region 3.

Further note that the segment boundaries are not smooth because of the
boundary being formed by the outermost Voronoi cells. The photons that
comprise the boundary are also subject to stochasticity, due both to PSF-
induced statistical variations in photon arrival locations, as well as the
greedy merging process. Visual inspection of the results suggests that at the
lowest surface brightness levels, fluctuations in the counts could result in
oversegmentation of what is usually considered the background (e.g., the two
large extended regions along the left side of the bottom edge of the
$\mathcal{F}$). Nevertheless, the expanded views of the inset regions in
Figure 7 show that the segmentation correctly separates diffuse emission
structures at different spatial scales and surface brightness levels. In
particular, transitions in the spatial density of photons across the
boundaries are clearly discernible by eye, such as those between segments
B$\leftrightarrow$C, B$\leftrightarrow$D, D$\leftrightarrow$E,
D$\leftrightarrow$F, E$\leftrightarrow$F in Region 1; between segment B and
segments A, C, D, F, G, I, J in Region 2; and segment A and segments B, C, D,
G, and H as well as C$\leftrightarrow$F and C$\leftrightarrow$D in Region 3.
Some transitions are too subtle to be visually recognizable (e.g.,
A$\leftrightarrow$B in Region 1, B$\leftrightarrow$E and B$\leftrightarrow$H
in Region 2, and D$\leftrightarrow$E in Region 3) but are required due to the
computed contrasts in the counts per unit area. Conversely, the brightness
transitions across B$\leftrightarrow$C$\leftrightarrow$D in Region 1,
B$\leftrightarrow$C$\leftrightarrow$K and
B$\leftrightarrow$I$\leftrightarrow$L in Region 2, and
A$\leftrightarrow$B$\leftrightarrow$C, F$\leftrightarrow$C$\leftrightarrow$D,
and D$\leftrightarrow$A$\leftrightarrow$G,H are apt demonstrations of the
capability of SRGonG to perform at the level of human visual acuity.
Parametric modeling to capture the spatial variations in such structures would
be much more difficult than the segmentations achieved here.

Figure 8: Segmentation of the photon list shown in Figure 5. Photons are
color-coded by the estimated intensity of their segment. Left Panel: 55 fitted
segments obtained with SRGonG with seeds specified via a grid supplemented by
local maxima. Right Panel: 61 fitted segments obtained for brute-force
segmentation, where all photons are assumed to be seeds. The differences
between these two segmentations are described in Section 6.1.

Figure 9: Identifying point-like sources in the event list data in Figure 5.
Left Panel: ”Point-like” SRGonG segments selected to have segment area
$\leq$0.0003 (31 total segments; 22 after concentric segments are merged).
Right Panel: Sources detected with wavdetect (Zezas et al., 2002) overlaid on
a full-band image of the Chandra data (OBSID 315). The size of the field is
$2.8\times 2.3$ arcmin. Red circles mark detected sources that were shown to
be inconsistent with Chandra’s PSF according to the analysis of (Zezas et al.,
2002), while cyan circles mark the remaining sources (some of which may also
be extended). The size of the circles corresponds to the optimal wavdetect
scale. Thicker solid circles indicate sources that are also identified by
SRGonG; the remaining sources are indicated by the thinner dashed circles. The
cyan and orange numbers give the source counts measured by the customized
aperture photometry of Zezas et al. (2002) and the SRGonG method respectively
(the SRGonG-based counts are not background subtracted). Note that most strong
point-like sources are identified as such by the SRGonG method, and in general
there is good agreement in the estimated source intensities.

An important factor in obtaining a reliable segmentation is the initial seed
specification. It is worthwhile establishing that the scheme we propose
generates a useful segmentation and does not miss features. For this, we
compare the SRGonG method against a brute-force segmentation where every
photon in the dataset is taken to be a seed, and the corresponding Voronoi
cells are merged using the BIC criteria as described in Section 3.3. This
brute-force scheme is similar to Scargle’s (2002) method (but using the BIC
criteria instead of Bayes Factors) in that it eschews the Seeded Region
Growing on Graph step developed and described in Section 3.1.

In Figure 8, we compare the SRGonG-based segmentation (left panel) against the
brute-force segmentation (right panel). Although at first glance the two
segmentations look similar, a closer inspection reveals crucial differences.
While the quality of the identification of point-like sources does not differ
significantly, there are significant differences in the diffuse emission
regions that strongly favor the SRGonG segmentation. Notice that segment C of
Region 1 from SRGonG (top panel of Figure 7) is missing in the brute-force
segmentation, and is effectively subsumed into segment D, which in turn also
subsumes segment B. These changes are prima facie unsupported by the visible
variations in the surface density of photons. Similarly, we see that segment E
is incorrectly extended, and a different segment extends down into the middle
of segment D. Such cases are also seen in Region 2 (middle panel of Figure 7),
where all of the complexity found as segments C, D, and E are lost in the
brute-force segmentation; and in Region 3 (bottom panel of Figure 7) where the
clear separation of segments A and B is lost in the brute-force segmentation,
as is the point-like source at $\sim$(0.4,0.15). In summary, clear variations
in surface brightness are recovered in the SRGonG segmentation, unlike in the
brute-force method. Using a smaller set of seeds improves the robustness of
the segmentation by avoiding the chaotic development of early merging steps;
errors in early stages accumulate because of the greedy merging process. We
thus conclude that SRGonG is superior because of, and not despite, the much
smaller, but perceptively selected number of seeds used to carry out the
segmentation.888Just as Markov Chain Monte Carlo techniques rely on running
multiple chains and verifying consistency to gain confidence in the analysis
results (Gelman & Rubin, 1992), we recommend that analyses that use SRGonG
also consider the sensitivity of the results to the adopted seed set. The
schemes that we recommend in Section 3.2 are adequate to handle most scenarios
encountered in astronomy, but are nonetheless characterized by several run-
time specified parameters ($m_{\rm grid}$, $m_{\rm graph}$, $m_{\rm nn}$,
$m_{\rm R}$, $m_{\rm vorthr}$). Work to formalize this process via bootstrap
analysis is ongoing (Jue Wu, private communication).

Although SRGonG is not designed as a point-source detection method, it is
instructive to see how it behaves in the case of point-like sources. In Figure
9 we how point-like sources can be identified in SRGonG (left panel) compared
to a wavelet-based method that is optimized to detect point sources (right
panel; wavdetect; Freeman et al., 2002). Based on the typical size of the
Chandra point spread function (PSF), we isolate all SRGonG segments that cover
an area of comparable or smaller size to the PSF999We choose segments
identified by areas ${{\rm Area}(\mathcal{S}_{k})}\leq{0.0003}$ in normalized
coordinates, which corresponds to a circular area of radius $\approx$1.6′′ on
the sky, comparable to the extent of the Chandra PSF.  and show them in the
left panel. (We emphasize that this is not a method to detect point-like
sources; while SRGonG segments with larger areas than the PSF size can be
flagged as extended, regions with small areas cannot be definitively flagged
as point sources, since such segments can occur due to layered segmentation of
extended sources or even due to statistical fluctuations in the surface
brightness of diffuse emission.) In the right panel, we show all the wavdetect
detected sources, superposed on a counts image of the same field. Since
wavdetect is optimized to find point sources in a variety of scales, it may
also detect more diffuse sources. Sources that are identified as extended
based on visual inspection and/or comparison with the PSF profile (Zezas et
al., 2002, e.g. lack of a core, or PSF fitting for sources with more than 100
counts;) are marked by red circles, while point-like sources are marked by
cyan circles. We note that the SRGonG segmentation invariably finds as “point-
like” (based on the segment area criterion) the point sources that are
confirmed by the inspection process of Zezas et al. (2002) and does not find
the extended sources identified by wavdetect. The latter are instead
components of larger diffuse emission segments. In this respect, although the
SRGonG is not a point-source detection algorithm, screening of the identified
segments based on the segment area, ${{\rm Area}(\mathcal{S}_{k})}$, can be
used to distinguish extended regions from point-like sources.

###  6.2 Advantages and Limitations of SRGonG

Our simulations (Section 4) and analysis of the Chandra Antennae dataset
(Sections 5 and 6.1) illustrate SRGonG’s strength in identifying sources at
many different scales. The method allows the identification of extended
diffuse structures in X-ray data regardless of their shape, i.e., no
assumptions are made about the morphology or the homogeneity of the sources.
Note that while the blurring due to the shape of the PSF is not explicitly
modeled, this has negligible effect on any source structure at scales larger
than the size scale of the PSF. Thus, we expect that useful results can be
obtained even when the PSF varies across the field of view, which can happen
due to several reasons (the quality of the telescope optics can degrade away
from the aimpoint; or fields are observed which have a large diversity of soft
and hard sources, each with significantly different PSF shapes and sizes; or
when complex combinations of datasets, such as multiple observations carried
out at different angular offsets, are combined). At spatial scales larger than
the PSF size, we expect that results are not reliant on the specific
characteristics of the PSF101010Note that PSF size information may be
incorporated into the analysis, e.g., by requiring that any segment that is
found to have a smaller area than that of the PSF be subsumed into a
surrounding or adjacent segment. We do not use such a criterion in this work,
though such a strategy is demonstrated in Figure 9..

Even when the photons are sparsely distributed, e.g., when the observation is
dominated by diffuse structures at low surface brightness, and blurring due to
the shape of the PSF is not included, point sources that may exist in the
field of view can be identified due to the increased concentration of photons
at their locations. However, note that SRGonG is designed to identify large
scale extended regions, so the focus and trade-offs are different. Thus, weak
point sources with low contrast against the surrounding diffuse emission are
likely to be subsumed into the diffuse regions. But because these are by
definition weak, they are unlikely to contribute significantly to the
brightness (or hardness) of the diffuse component. This situation is
effectively similar to the situation where the detection sensitivity of a
telescope is insufficient to resolve apparently diffuse emission into its
point source population. An additional issue to consider is the bias in the
point source brightness demonstrated in Figure 3. This bias arises as area
fluctuations from small point sources are naturally bounded at zero, but can
extend into the area of the diffuse emission, leading to a skewed error
distribution. So point source intensities found by SRGonG should not be used
directly, but must be re-estimated using appropriate techniques (e.g., Primini
& Kashyap, 2014). However, note that the bias demonstrated in Figure 3 comes
from point sources that do not have PSF wings; in real sources where point
sources are sharply peaked due to the PSF, the area distribution bias works to
the advantage of SRGonG, incorporating more of the PSF wings into the point
source and reducing the resulting contamination of the diffuse emission by
strong point sources.

Unlike adaptive smoothing, source detection, or contouring methods, SRGonG
does not set S/N thresholds or rely on thresholds of source significance to
determine the presence or extent of contiguous regions. Thus, even regions
that may be characterized by low surface brightness tend not to be over-
segmented. Conversely, since the uncertainty in the estimated brightness is
dependent on the number of photons that fall within the segment, small
variations in adjacent regions can be more easily distinguished when the areas
of the segments are sufficiently large.

Also of note is that SRGonG works directly on photon lists, the most basic
form of high-energy X-ray and $\gamma$-ray datasets and the Poisson nature of
the data is explicitly accounted for during the merging process. While this
can have detrimental effects on running time when the size of the dataset is
large111111For illustration, the analysis of the Antennae dataset, with
$50,700$ photons and $491$ seeds, takes $\approx$240 s on a 2021 epoch 14′′
MacBook Pro with an Apple Silicon M1 SOC., using the data at the highest
available resolution avoids the requirements to define artificial binning
sizes.

Of greater concern is the dependence of SRGonG results on the distribution of
the initial seeds, especially for fields with low contrast. This may result in
unstable behaviour because of fluctuations in the local minima in the spatial
intensity distribution, leading to both false segmentation and false merging.
We caution that while the schemes we describe in Section 3.2 are generally
adequate and perform well (see, e.g., Section 6.1), as is typical with seeded-
region-growing methods, the sensitivity of the segmentation to the adopted
seed structure must always be checked.

Future extensions of this method will include quantification of the
uncertainty of the segmentation resulting from the stochastic nature of the
data, which would be quantified in terms of uncertainty on the number of
segments, the outline of the segments, and the corresponding source flux
within each segment. Other avenues to explore include different merging
procedures as substitutes for the greedy merge to address the over-
segmentation. The goal of such alternative merging options would be to search
more possible final segmentations and make the algorithm more robust to seed
initialization. Yet another potential extension is to perform the analysis in
3 dimensions, incorporating photon energy information. Currently spectral
information can only be used by running the code on passband filtered data.

##  7 Summary

We have developed an algorithm that provides a piece-wise constant
segmentation of a photon event list that approximates the spatial structure
present in the data. Point-wise surface brightnesses are initially estimated
as the inverse of their Voronoi cell areas and cells with similar brightnesses
are grouped together to grow segments. The seeds needed to grow the segments
can be initialized as regular grids, additionally supplemented with local
maxima, or set using more complex processes by stratified sampling of Voronoi
cell areas. The process begins with a deliberate over-segmentation, and
neighboring segments are sequentially merged by maximizing the BIC change. The
resulting (greedy) segmentation generates apertures on the sky plane that can
be used to collect photons and carry out further analysis in a way that
removes manual intervention in selecting regions of interest. We have explored
this method via both simulations and application to a complex Chandra dataset,
and find that it consistently provides a good description of both point-like
and extended diffuse regions of arbitrary shapes.

We note that this is not a source detection method, but a robust method for
the definition of source regions, especially for extended sources. In this
way, it can be used to perform photometry or spectroscopy on arbitrarily
shaped extended sources.

This method provides several advantages over other commonly used methods for
the analysis of extended sources in high-energy photon data. Namely, it allows
the identification of sources at different scales even when they are embedded
within each other without imposing any restrictive assumptions on the spatial
distribution of the source photons or the source intensity.

(ACIS)

(https://cxc.harvard.edu/ciao/; Fruscione et al., 2006); PINTofALE (Kashyap &
Drake, 2000); Matlab (https://www.mathworks.com/products/matlab.html); SRGonG
(https://github.com/jujWang96/Astro_sim).

## Acknowledgements

This work was conducted under the auspices of the CHASC International
Astrostatistics Center. CHASC is supported by NSF DMS-18-11308, DMS-18-11083,
DMS-18-11661, DMS-21-13615, DMS-21-13397, and DMS-21-13605; by the UK
Engineering and Physical Sciences Research Council [EP/W015080/1]; and by NASA
18-APRA18-0019. We thank our CHASC colleagues for many helpful discussions,
especially Jilei Yang for his valuable comments on an earlier draft. MF, JW
and TCML acknowledge further support from NSF through CCF-19-34568,
DMS-18-11405 and DMS-19-16125. VLK further acknowledges support from NASA
contract to the Chandra X-ray Center NAS8-03060. DvD and AZ were also
supported in part by a Marie-Skodowska-Curie RISE (H2020-MSCA-
RISE-2015-691164, H2020-MSCA-RISE-2019-873089) Grants provided by the European
Commission.

## References

  * Adams & Bischof (1994) Adams, R., & Bischof, L. 1994, IEEE Transactions on Pattern Analysis and Machine Intelligence, 16, 641
  * Akhlaghi & Ichikawa (2015) Akhlaghi, M., & Ichikawa, T. 2015, ApJS, 220, 1
  * Allard & Fraley (1997) Allard, D., & Fraley, C. 1997, Journal of the American Statistical Association, 92, 1485
  * Angel et al. (2010) Angel, S., Parent, J., & Civco, D. L. 2010, The Canadian Geographer/Le Géographe canadien, 54, 441
  * Aue & Lee (2011) Aue, A., & Lee, T. C. M. 2011, The Annals of Statistics, 2912
  * Baldi et al. (2006) Baldi, A., Raymond, J. C., Fabbiano, G., et al. 2006, ApJS, 162, 113
  * Barr & Schoenberg (2010) Barr, C. D., & Schoenberg, F. P. 2010, Biometrika, 977
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, Astronomy & Astrophysics Supplement, 117, 393
  * Bobin et al. (2016) Bobin, J., Sureau, F., & Starck, J.-L. 2016, A&A, 591, A50
  * Chiu et al. (2013) Chiu, S. N., Stoyan, D., Kendall, W. S., & Mecke, J. 2013, Stochastic geometry and its applications (John Wiley & Sons)
  * Ebeling et al. (2006) Ebeling, H., White, D. A., & Rangarajan, F. V. N. 2006, Monthly Notices of the Royal Astronomical Society, 368, 65
  * Ebeling & Wiedenmann (1993) Ebeling, H., & Wiedenmann, G. 1993, Physical Review E, 47, 704
  * Fabbiano et al. (2001) Fabbiano, G., Zezas, A., & Murray, S. S. 2001, ApJ, 554, 1035
  * Fan & Lee (2014) Fan, M., & Lee, T. C. M. 2014, IET Image Processing, 9, 478
  * Farias et al. (2020) Farias, H., Ortiz, D., Damke, G., Jaque Arancibia, M., & Solar, M. 2020, Astronomy and Computing, 33, 100420.  https://www.sciencedirect.com/science/article/pii/S2213133720300743
  * Freeman et al. (2002) Freeman, P. E., Kashyap, V., Rosner, R., & Lamb, D. Q. 2002, The Astrophysical Journal Supplement Series, 138, 185
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Observatory Operations: Strategies, Processes, and Systems, Vol. 6270, International Society for Optics and Photonics, 62701V
  * Gelman & Rubin (1992) Gelman, A., & Rubin, D. B. 1992, Statistical Science, 7, 457 .  https://doi.org/10.1214/ss/1177011136
  * Gondhalekar et al. (2022) Gondhalekar, Y., de Souza, R. S., & Chies-Santos, A. L. 2022, Research Notes of the American Astronomical Society, 6, 128
  * Hausen & Robertson (2020) Hausen, R., & Robertson, B. E. 2020, ApJS, 248, 20
  * Hubert & Arabie (1985) Hubert, L., & Arabie, P. 1985, Journal of classification, 2, 193
  * Jones et al. (2015) Jones, D. E., Kashyap, V. L., & Van Dyk, D. A. 2015, The Astrophysical Journal, 808, 137
  * Kashyap & Drake (2000) Kashyap, V., & Drake, J. J. 2000, Bulletin of the Astronomical Society of India, 28, 475
  * Lee (2000) Lee, T. C. M. 2000, Journal of the American Statistical Association, 95, 259
  * Magnussen et al. (2006) Magnussen, S., Allard, D., & Wulder, M. A. 2006, Scandinavian Journal of Forest Research, 21, 239
  * Meyer et al. (2021) Meyer, A. D., van¬†Dyk, D. A., Kashyap, V. L., et al. 2021, Monthly Notices of the Royal Astronomical Society, 506, 6160.  https://doi.org/10.1093/mnras/stab1456
  * Møller (1994) Møller, J. 1994, Lectures on random Voronoi tessellations, Vol. 87 (Springer Science & Business Media)
  * Otsu (1979) Otsu, N. 1979, IEEE transactions on systems, man, and cybernetics, 9, 62
  * Peng et al. (2011) Peng, B., Zhang, L., & Zhang, D. 2011, IEEE Transactions on Image Processing, 20, 3592
  * Picquenot et al. (2019) Picquenot, A., Acero, F., Bobin, J., et al. 2019, A&A, 627, A139
  * Picquenot et al. (2021) Picquenot, A., Acero, F., Holland-Ashford, T., Lopez, L. A., & Bobin, J. 2021, A&A, 646, A82
  * Primini & Kashyap (2014) Primini, F. A., & Kashyap, V. L. 2014, ApJ, 796, 24
  * Rand (1971) Rand, W. M. 1971, Journal of the American Statistical Association, 66, 846.  https://www.tandfonline.com/doi/abs/10.1080/01621459.1971.10482356
  * Sanders (2006) Sanders, J. S. 2006, MNRAS, 371, 829
  * Sanders & Fabian (2001) Sanders, J. S., & Fabian, A. C. 2001, MNRAS, 325, 178
  * Scargle (2002) Scargle, J. D. 2002, in AIP Conference Proceedings, Vol. 617, AIP, 163–173
  * Scargle et al. (2013) Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, The Astrophysical Journal, 764, 167
  * Starck et al. (2002) Starck, J. L., Pantin, E., & Murtagh, F. 2002, PASP, 114, 1051
  * Vikhlinin et al. (1998) Vikhlinin, A., McNamara, B. R., Forman, W., et al. 1998, ApJ, 502, 558
  * Zezas et al. (2006) Zezas, A., Fabbiano, G., Baldi, A., et al. 2006, ApJS, 166, 211
  * Zezas et al. (2002) Zezas, A., Fabbiano, G., Rots, A. H., & Murray, S. S. 2002, ApJS, 142, 239

##  Appendix A Nearest Neighbor Labeling

Here we describe the heuristic by which selected photons are collected into
groups characterized by their proximity (used in the Voronoi-area stratified
sampling scheme for seed specification; see Section 3.2). The photons
considered in a given stratum are defined by a small range of Voronoi areas,
or analogously, are located at similar contour levels if an image were
constructed from the photons. Thus, they are likely to be sparsely
distributed, but with clumps of photons surrounding higher intensity regions.
The goal here is to group the clumped photons that are near each other,
without breaking up rings or other complex shapes. We emphasize that this
heuristic is a quick but approximate pre-processing method to pick seeds for
the full-fledged SRGonG algorithm. We expect this heuristic to be useful in
situations where the astronomical dataset is characterized by sparsely
distributed structures with a large dynamic range in surface brightness.

We first determine an average characteristic length scale for the ensemble of
photons included in stratum $\Upsilon$, as

${{L}}_{\Upsilon}=2\sqrt{\frac{1}{2}\left[\max_{i\in\Upsilon}\\{{{\rm
Area}(\mathcal{V}_{i})}\\}+\min_{i\in\Upsilon}\\{{{\rm
Area}(\mathcal{V}_{i})}\\}\right]}.$

This ensures that the length scale is typical of stratum $\Upsilon$. We begin
with an arbitrary photon from $\Upsilon$, assigning it a unique group label,
and recursively assign this group label to any neighbor, i.e., any photon in
stratum $\Upsilon$ located within a Euclidian distance of $L_{\Upsilon}$ from
any photon assigned to this group. The recursive labeling ends when no new
neighbors are present, and we move to another arbitrary as yet unlabeled
photon in $\Upsilon$, assign it a different label, and repeat the process. We
continue this labeling until all photons in $\Upsilon$ are assigned labels.
For a case where the photons are placed uniformly on a regular grid, this
results in all the photons being aggregated into one clump with one label. If
there are multiple clumps separated by $>L_{\Upsilon}$, each clump will be
assigned a separate label. We eventually discard all clumps with fewer than
$m_{\rm vorthr}$ photons and do not use them to set a seed. The entire process
is repeated for each of the $m_{\rm R}$ strata.


# Localization using Mobile Wireless Sensor Networks

Aaron John-Sabu  Dept. of Electrical Engineering
Indian Institute of Technology Bombay
Mumbai, India
<EMAIL_ADDRESS>
###### Abstract

Wireless Sensor Networks (WSNs) are groups of spatially distributed and
dedicated autonomous sensors for monitoring (and recording) the physical
conditions of the environment (and organizing the collected data at a central
location). They have been a topic of interest due to their versatility and
diverse capabilities despite having simple sensors measuring local quantities
such as temperature, pH, or pressure. We delve into understanding how such
networks can be utilized for localization, and propose a technique for
improving conditions of living for animals and humans on the IIT Bombay
campus.

###### Index Terms:

Mobile Wireless Sensor Networks, localization, target tracking, sustainable
living

##  I INTRODUCTION

In a society that is booming with collaborative and wide-spread activities,
both human-based and autonomous, it is a vital requirement to incorporate
similar networks that can monitor and support human needs such as pollution
control, navigation, industrial monitoring etc. Utilizing wireless sensors
networks for such functionalities further helps in the reduction of power
consumption, work time, and difficulties in implementation. We study the work
performed on the application of WSNs in localization. The learnings are
conglomerated to develop a mechanism for the maintenance of wildlife in the
IIT Bombay campus, particularly dogs.

##  II PREVIOUS RESEARCH

In recent times, the incorporation of WSNs into populated societies have
improved conditions of life and work, from increasing crop yield and
optimizing utility distribution systems (electricity, water, etc.) to reducing
wasted power and avoiding traffic jams. [1] elucidates how WSNs have developed
throughout the years, from more expensive and fewer sensors to cheaper but
larger quantities of sensors, from discrete circuits and multi-chip solutions
to system-on-chip (SoC) devices, and from one-way communication links to
bidirectional links and mesh and star designs. These optimizations have led to
widespread applications of WSNs one of them being navigational localization of
objects and bodies using methods such as trilateration.

###  II-A Localization using Wireless Sensor Networks

Localization can be performed using multiple methods, a very common one being
GPS. However, the accuracy of GPS data is relatively low ($>10m$). [2]
proposes a more accurate localization method wherein a wireless sensor network
of ZigBees is developed and their relative signal strengths are used for
trilateration-based localization. The Mobile Target Tracking (MTT) problem
intends to find the moving path of a target in a field based on target
locations that are sampled at random intervals. [3] discusses several
algorithms for solving this problem using two aspects: determining the current
location of the target (localization, path tracing), and processing
information collaboratively among multiple sensor nodes. Traditional methods
involving the informed selection of sensors, binary sensor-based methods with
centralized and distributed architectures, and other methods based on
triangulation are suggested for tracking, while information can be processed
using leader-based algorithms or distributed algorithms.
Unlike open environments, locations with several obstructions or jamming
hinder the proper function of the Global Position System (GPS). In such
scenarios, it is necessary to develop a positioning system that can complement
GPS. [4] proposes a pedestrian navigation system (PNS) wherein heterogeneous
agents are used for sensing, communication and relaying information. Here,
sensors such as accelerometers and magnetometers are used as part of a dead
reckoning (odometry) mechanism to localize the user. The NavMote (sensor on
the user) exchange information with NetMotes (predetermined sensors) when they
are in range; otherwise it works on its own. This system has been tested in
both indoor and outdoor environments, giving satisfactory performance at a
distance accuracy of $\pm 1\%$ and a heading accuracy of $1^{\circ}$.
Localization for mobile targets can be performed using beacon-based methods,
wherein some beacons aware of their positions provide geographic information
to ordinary sensor nodes to localize, whose precision increases with the
number of beacons. [5] proposes an algorithm utilizing mobile beacons that
traverse the network deployment area and broadcast beacon packets to generate
a number of virtual beacons. The distance between the sensor node and the
beacon can be calculated using RSSI (Received Signal Strength Indicator).
Since a single-mobile-beacon system can have issues such as co-linearity due
to the straight line moving trajectory of the mobile beacon, a three-mobile-
beacon-assisted mechanism is suggested wherein the object (sensor node,
$S_{i}$) is localized as the weighted centroid of the three beacon positions:

$S_{i}(\hat{x}_{si},\hat{y}_{si})=\frac{\sum_{j=1}^{m}w_{ij}V_{j}(x_{vj},y_{vj})}{\sum_{j=1}^{m}w_{ij}},w_{ij}=\frac{1}{(d_{ij})^{g}},$

where $V_{j}$ is the $j$th virtual beacon, $d_{ij}$ is the distance between
$S_{i}$ and $V_{j}$, and $g$ is an adjustable degree.
While anchors in anchor node-based schemes acquire their positions in advance
using GPS systems or artificial arrangement to locate unknown nodes, unknown
nodes in anchor node-free schemes are located using the connectivity
information between unknown nodes and anchor nodes. The former achieves better
localization accuracy but require a large number of anchor nodes which
increases the energy consumption and hardware cost of WSNs. [6] proposes using
mobile anchor nodes to maximize the localization accuracy while decreasing the
energy consumption of WSNs. An anchor node moves based on an equilateral
triangle trajectory in a WSN area and broadcasts position and time messages
periodically which, on reception, are used via RSSI-based trilateration to
determine the position of unknown nodes. Though sensitive to the standard
deviation of noise, the algorithm reduces the number of beacon positions,
trajectory lengths and node density. Moreover, it remains robust at high
traveling speeds of the anchor node.

###  II-B Optimization of WSNs and Other Applications

The distributed detection process involves sensor nodes that are deployed
randomly in a field to collect data from the surrounding environment. [7]
describes and compares three detection schemes: centralized, distributed, and
quantized, based on parameters such as detection probability and overall
probability of error, measured against varying energy consumption, time
periods and distances from the sensors to the command center. The robustness
of each scheme is also tested using node destruction and partial observation
deletion. The distributed scheme is observed to be superior in energy
consumption (especially over large distances) and robustness. Although the
centralized scheme uses fewer nodes, the distributed option needs fewer than
twice that number to achieve the same detection performance.
The neighbor discovery problem deals with situations where sensor nodes find
their neighbors constantly in mobile sensor networks by communicating with
each other while in motion and forward the collected information to a central
command center. Here, the active and dormant status of sensors can be
controlled, hence reducing the energy consumption significantly. However, this
may cause additional discovery latency since discovery is possible only when
neighboring nodes have overlapping active slots. Previous research had
introduced the group-based method where a third state for waking up actively
is used to communicate the schedule and verify the neighborhood of nodes. [8]
proposes an algorithm that considers the embedded spatial properties and
actively modifies the active time of nodes depending on the number of
undiscovered neighbors. This has been tested using simulations and the
discovery time has been found to be minimal when compared to algorithms
presented in existing literature.
[9] proposes a wireless sensor network distributed in the physical environment
for sensing, actuation and communication such that indoor lighting is
automated, while also giving override capabilities to users based on
individual preferences. The system is decomposed into smaller pseudo-static
subnets and issues such as nodes belonging to multiple subnets are solved
using a multi-agent system approach wherein agent actions are averaged or
voting can be used to determine the action. Furthermore, a supervised learning
model can be used to predict the action of the user and act accordingly while
remaining bounded to physical constraints. This can be replaced by model-free
reinforcement learning algorithms so that the agents are trained when a model
is unavailable or too difficult to compute in a complicated real-world
environment.

##  III The Effect of TAG LINKAGES for MOBILE WSN-based LOCALIZATION

###  III-A Introduction

The Indian Institute of Technology Bombay campus, spread across over 500
acres, is home not only to humans but also a wide array of plants and animals
from leopards and crocodiles to cows, dogs and cats. Unfortunately, in recent
times, human activities in the institute have been disrupted by stray dogs.
Also, survival for these dogs has also become difficult due to the movement of
heavy-duty construction vehicles along the same routes posing a threat to
their lives. This report proposes a technique to track and guide dogs while
avoiding harm to the human population as well as the dogs in the institute.

###  III-B Procedure and Solution

Based on the concept of mobile target tracking as discussed in [3], a network
of distributed sensors may be placed at suitable locations in the institute.
However, owing to the large size of the campus, this is not a scalable
solution since the number of sensors and the corresponding energy consumption
will be huge. Moreover, the distribution of dogs across the institute need not
be uniform, due to which placing sensors at certain locations will not be
efficient although there is a possibility for a small number of dogs to visit
these areas.
Incorporating ideas from [6], the movable sensors on the dogs can behave as
anchors using which other movable sensors can localize themselves via
trilateration. The number of fixed sensors can hence be reduced if the number
of dogs is large. This is demonstrated using a simulation wherein the number
of fixed sensors has been reduced to $6$ in a $10m\times 10m$ area where each
sensor can sense sensors within a radius of $7m$. Three movable sensors are
placed such that the size of their point represents the uncertainty of finding
their location.
We define the neighbors of a sensor at a particular moment to be the sensors
capable of connecting to and providing information about position and time to
the sensor for trilateration at that instant. A simple geometric uncertainty
principle is used wherein the uncertainty of the location of the sensor
increases given the number of neighbors is less than $4$ (with a lower rate
for $3$ neighbors) and decreases otherwise at an exponential rate proportional
to the number of neighbors.
Free-ranging dogs generally exhibit territoriality as studied in [10] and
those in the institute have been observed to do so in considerably small
pieces of land such as individual hostels and small strips of roads.
Similarly, the three sensors in the simulation are constrained within
particular regions.

###  III-C Results

The mechanism is simulated in two-dimensional space within the given
dimensions. We initialize the sensors with some positional uncertainty as
shown in Fig. 1.

Figure 1: Initial Configuration - Without Tag Links

On deactivating links between tags (i.e., they do not behave as mutual
anchors), there are larger periods when the number of links is $2$ due to
which the uncertainty of some sensors is observed to diverge with time as
depicted in Fig. 2.

Figure 2: Final Configuration - Without Tag Links

The second simulation involves links between neighboring tags. This starts
from the beginning of the simulation as shown in Fig. 3.

Figure 3: Initial Configuration - With Tag Links

For simplicity, we assume that the range of the tag sensors is equal to that
of the anchors although, in a practical scenario, this might not be the case
since the former will consist of simpler circuits and transceivers, hence
resulting in smaller ranges. It is observed from Fig. 4 that the uncertainty
of the previously mentioned sensors does not diverge since the sensors are in
contact with each other for several periods which previously had led to the
divergence of uncertainty.

Figure 4: Final Configuration - With Tag Links

##  IV FUTURE SCOPE

Wearing traditional sensors for a long duration may be detrimental to the
health of the dog, and the sensing quality may degrade from the environmental
and hygienic conditions. Biosensors and flexible electronics may be more
appropriate long-term options. The long-term effects of biologic cybernetics
and electronic stimulation on dogs is open to research although there has been
progress in short-duration studies on several animal species as shown in [11],
[12] and [13].

The simulation code related to this project is available at
github.com:aaronjohnsabu1999/Localization-Using-WSNs

##  V ACKNOWLEDGMENTS

I am grateful to Prof. Siddharth Tallur for the guidance and support
throughout this project as the course instructor for EE 617 (Sensors in
Instrumentation) at IIT Bombay.

## References

  * [1] A. Benefit, “The evolution of wireless sensor networks”, [online] Available: https://www.silabs.com/documents/public/white-papers/evolution-of-wireless-sensor-networks.pdf
  * [2] Lee, J-H, Kim, K, Lee, S-C, and Shin, B-S, “An Efficient Localization Method Based on Adaptive Optimal Sensor Placement,” Int. Journal of Distributed Sensor Networks, 2014
  * [3] Gupta, A, Gui, C, and Mohapatra, P, “Mobile Target Tracking using Sensor Networks,” Mobile, Wireless, and Sensor Networks: Technology, Applications, and Future Directions, 2006
  * [4] Fang, L, Antsaklis, PJ, et al., “Design of a Wireless Assisted Pedestrian Dead Reckoning System—The NavMote Experience,” IEEE Transactions on Instrumentation and Measurement, 2016
  * [5] Cui, H-q, Wang, Y-l, et al, “Three-Mobile-Beacon Assisted Weighted Centroid Localization Method in Wireless Sensor Networks,” Int. Conf. on Software Engineering and Service Science, 2011
  * [6] Jiang, J, Han, G, Xu, et al., “LMAT: Localization with a Mobile Anchor Node Based on Trilateration in Wireless Sensor Networks,” IEEE Global Telecommunications Conf., 2011
  * [7] Yu, L and Ephremides, A, “Detection, Energy, and Robustness in Wireless Sensor Networks,” Mobile, Wireless, and Sensor Networks: Technology, Applications, and Future Directions, 2006
  * [8] Niu, Q, Bao, W, and Xia, S, “An Improved Group-Based Neighbor Discovery Algorithm for Mobile Sensor Networks,” Int. Journal of Distributed Sensor Networks, 2014
  * [9] Sandhu, JS, Agogino, AM, and Agogino AK, “Wireless sensor networks for commercial lighting control: Decision making with multi-agent systems,” AAAI Workshop on Sensor Networks, 2004
  * [10] Majumder, SS, et al., “A dog’s day with humans – time activity budget of free-ranging dogs in India,” Current Sciences 10(6), 2014
  * [11] Johnson, LA, and Fuglevand, AJ, “Mimicking muscle activity with electrical stimulation,” Journal of Neural Engineering vol. 8(1), 2011
  * [12] Rezaee, Z, and Kobravi, HR, “Human Gait Control Using Functional Electrical Stimulation Based on Controlling the Shank Dynamics,” Basic and Clinical Neuroscience vol. 11(1), 2020
  * [13] Cao, F, Zhang, C, Vo Doan, TT, et al, “A Biological Micro Actuator: Graded and Closed-loop Control of Insect Leg Motion by Electrical Stimulation of Muscles,” PLOS ONE 9(8), 2014