Visualizing the Wavenumber Content of a Point Pattern

Spatial point patterns are a commonly recorded form of data in ecology, medicine, astronomy, criminology, epidemiology and many other application fields. One way to understand their second order dependence structure is via their spectral density function. However, unlike time series analysis, for point patterns such approaches are currently underutilized. In part, this is because the interpretation of the spectral representation of the underlying point processes is challenging. In this letter, we demonstrate how to band-pass filter point patterns, thus enabling us to explore the spectral representation of point patterns in space by isolating the signal corresponding to certain sets of wavenumbers.


Introduction
Spectral analysis of time series and random fields is a very mature research area [1,2].Spatial point patterns are found in a wide ranging set of applications, see for examples [3,4,5].In the point process setting, some work has been done to establish spectral representations [6] and estimation of the spectral density function [7,8,9], but the methodology has yet to be widely adopted.In part, this is because interpreting the spectral representation of point processes is much harder than interpreting the spectral representation of random fields or time series, a problem we seek to address in this paper.
One way to think about the spectral density function is that it describes how the variability of the process can be explained by phenomena occurring at different wavenumbers.In order to better understand such phenomena, it can be useful to isolate the spatial variability associated with a certain selection of wavenumbers and relate this to the features of the point process that this selection represent.In a time series context, when these wavenumbers are an interval, this projection is called a band-pass, low-pass or high-pass filter, depending on the form the interval takes.Indeed, the earliest forms of spectral estimates were constructed by first band-pass filtering a process and then computing the variance of the output [1].
In our case, we wish to extend such an approach of wavenumber localization to spatial point processes.Linear filters for point processes have already been developed [6], and have been used for operations such as binning [10] or studied as an inverse problem, i.e., recover-original points ing points from a filter output [11].In contrast, we use such filtering to focus on the spatial behavior explained by different wavenumber ranges, a technique illustrated in Fig. 1.
In [9] we proposed efficient methods of spectral estimation, but the interpretation of a spectral decomposition needs to be established for point processes.This wavenumber interpretation of point processes is challenging, in part because points are essentially the roughest process one can construct, in other words, the furthest from a regular wave.In the wavenumber domain, this corresponds to having significant variability present at all wavenumbers, in contrast to random fields, which are often assumed to be band-limited.To address this, we show how to first project the observed point pattern into a suitable linear space, such that we retain the part of the signal associated with a given set of wavenumbers, and then, via illustrative examples, discuss how this adds insight to our understanding of such processes.

Background
We begin with some brief background on point processes.A point process X can be thought of as a random collection of locations in some space, say R d .A formal way to characterize such a process is by the function which counts the number of points in a given region B in space.This function is called the counting measure [6].A point process is said to be homogeneous if its stochastic properties are invariant to spatial shifting of the set B. As in the case of time series [12], homogeneous point processes have a spectral representation [6].Analogously to time series, the spectral representation of a point process is usually given with the mean removed.In particular, [6] define the mean corrected counting measure N 0 (B) = N(B) − λℓ(B), where λℓ(B) = E [N(B)], for any B.Here λ is called the intensity of the point process, and ℓ(•) is the Lebesgue measure.The spectral representation theorem [6] states for ϕ(•) belonging to a certain class of functions, and where We refer to Z(•) as the spectral process.
A consequence of homogeneity is that the spectral process is an orthogonal increment process, meaning that it is uncorrelated when evaluated over two disjoint sets.However, the variance of the spectral process is not constant: we have traded uncorrelatedness for inhomogeneity.Under certain conditions (see [6]), the point process in question has a spectral density function f (•), which satisfies In other words, the spectral density function describes the variance of the spectral process Z(•), and provides a decomposition of variance across scales.
In order to understand the spectral density function, it can therefore be helpful to understand the spectral process.Such a process cannot be visualized directly in a convenient fashion, as it is complex-valued and extremely noisy.However, we can make a projection of the point process of interest, so that we preserve the associated spectral process at a certain limited set of wavenumbers.This has the added benefit of producing a spatial domain process, which can be easier to interpret.In the language of signal processing, a subset of such projections would be called band-pass filters.Note that this is not a new way of estimating the spectral density (see our work [9]), but merely a way of separating the spectral process into different regions of the wavenumber space, and visualizing the result in the original space.

Band-Pass Filtering in Space
Let K ⊂ R d be a region in wavenumber space.Our task is to map the point process X to some new process Y(•) whose spectral process is the same as the spectral process of X in K, and zero outside.Let the filtered process of X be where h K (u) = K e 2πik•u dk.In the language of signal processing, h K (•) is the impulse response function.We show in Appendix A that Therefore, Y(•) is precisely the process we wanted.The mean term is present when the zero wavenumber is included in the region, but is easily removed if a mean zero process is desired.
Notice that the process Y(•) is not a point process, but a random field.It is possible to construct linear filters which map point processes to other point processes, see the example in [13].However, such processes cannot be concentrated on a bounded region in wavenumber, as point processes require variance from phenomena occurring at all wavenumbers [6].
In practice, we must choose the region K. Whilst there are multiple choices one might consider, we prefer annuli, as they preserve the directional properties of the underlying process (see Appendix B for more details).Typically, we are interested in behavior that is not completely random, in other words when the process does not behave like a Poisson process.Poisson processes have f (k) = λ for all k (they are the point process equivalent of white noise [6]).Therefore, we choose the specific wavenumber region by first estimating f (•), using the methodology of [9], and then choosing an annulus in which f (•) is substantially different from λ.

Applications
To illustrate our methodology, we explore both simulated and real data which are motivated by forest ecology.However, point patterns are an important form of data in many different applications.In particular, log-Gaussian Cox processes, a popular model in forest ecology [14], have also been utilized to model bacteria [15], seismic activity [16] and galaxies [17], among other applications.Furthermore, our technique can be applied to a wide range of other kinds of processes, such as Thomas processes [18], hyper-uniform processes [19] or any other homogeneous process, such as those found in [20].

Simulated Log-Gaussian Cox Processes with Noise
We begin by considering a simulated example in which we generate two point patterns with clustering occurring in the same locations, and then add observational noise to each in the form of Poisson processes, as shown in Fig. 2. Whilst it is clear from the clustered points that there is a joint behavior between the two processes, this is much less clear when we add the noise.However, by applying a low-pass filter, focusing on the wavenumbers at which the spectral density is substantially not λ, we can clearly see this correlation again.To generate this example, joint clustering is obtained by generating two log-Gaussian Cox processes [21] with shared intensity.Note that, whilst a single realization of such a process looks inhomogeneous, the random intensity is itself a homogeneous process, and thus the process as a whole is homogeneous.

Lansing Woods
In Fig. 3 we apply this filtering to a subset of the Lansing Woods data, recorded by [22].For the purpose of illustration, we focus on hickories and maples, which have been shown to have a spatially segregated relationship using established point process methods [14,23].We begin by computing spectral estimates from the point pat-terns (using the isotropic method proposed by [9]), and notice that the spectra are substantially larger than the processes' intensities for wavenumbers below a certain threshold.This indicates "non-Poisson" behavior associated with those wavenumbers.We therefore filter both processes with a low-pass filter, designed to isolate this behavior.We see in the filtered processes that this behavior corresponds to large scale clustering in both processes.Furthermore, these clusters are present at different locations in space, as expected and might indicate different soil type requirements and/or direct competition between individual trees of these species.

Computational Details
The band-pass filter proposed in Section 3 assumes that we observe the process on the entire space R d .Of course this is not possible in practice, where we will be limited to observations in some bounded region B. As a result, we must instead compute Therefore, we are missing contributions from points outside the observational window. 1 However, because the impulse response functions decay fairly quickly, Y(u) and Y B (u) are only different near the boundary of B.
Furthermore, the effect is often negligible relative to the other structures that are present.
Boundary issues withstanding, we can obtain very good properties in the wavenumber domain, obtaining near ideal band-pass filters in simulations.In the case of time series and random fields, it is not usually possible to apply ideal band-pass filters [1].However, in the special case of point processes we get close, as we can record the point process continuously (within a given region) without information loss, due to the sparse nature of points.This means that, up to some edge effects, we can digitally manipulate a continuous signal, resulting in sharp attenuation in wavenumber space.
In addition, our technique requires that the process is homogeneous.Furthermore, since it is merely a transformation of an observed realization, it is itself a stochastic process, and not an estimator converging to some true value.As such, one must be careful when interpreting the output.In particular, it is best used as a diagnostic and comparative tool, alongside other techniques such as spectral estimation [9], or to aid interpretation.

Conclusion
The analysis of point patterns differs from that of other stochastic processes, such as time series and random fields.One of these differences is that the wavenumber behavior of points needs to be present at all frequencies.To address this here, we have developed band-pass filtering for spatial point processes.We have shown how such an approach can be used to gain understanding of the spectral representation of a point process, exploring a variety of different models and processes of interest.We envision that this technique will provide greater interpretability for wavenumber domain analysis techniques, enabling them to be applied to a wide range of datasets.

Data Availability
An implementation of the methodology is available in the Julia package PointProcessFilters.jl [24], and the R package pppgram [25].Code to generate the figures in this paper are available on GitHub [26], and the Lansing Woods data is available from the R package spatstat [27].

A Filtering in the Wavenumber Domain
Let Ỹ(u) = Y(u) − λ1 K (0) be the mean corrected filtered process, then as required.2

B Choice of Region
The concept of a band in one-dimension has multiple generalizations in higher dimensions.However, the most obvious are annuli and shifted circles, though the box equivalents may be useful in some circumstances, see Fig. 4. We prefer the first of these shapes as it is invariant to rotation and does not influence the directional properties of the original data, indeed, standard Gaussian filters [28] used for band-pass filtering in image processing are also aiming to filter to an annulus.The construction of h K (•) for these regions is given in Appendix C. As discussed the choice of the specific annulus should be based on regions in which the spectral density function differs substantially from λ. Often this can be chosen by eye, but one may require formal tests in some settings (e.g.some modification of [29]), though this is beyond the scope of this paper.

C Constructing Impulse Responses
In order to construct the required filters, we need only the following relations.

C.1 Shifting
Consider applying a symmetric shift s to the region K, such that K = (K + s) ∪ (K − s) & (K + s) ∩ (K − s) = ∅ then from [30], Of course, one could shift without the symmetry, but the resulting process would be complex valued in general.

C.2 Unions
By additivity of integrals

C.3 Set Minus
Again, by additivity of integrals

Figure 1 :
Figure 1: Illustration of band-pass filtering of a point process.Here the points are a realization of a Thomas process, and the spectral density function (which is isotropic) is shown with a 1d slice.The filtered processes are plotted with darker shades indicating larger values, with a shared color range.

Figure 2 :
Figure 2: Illustration of the signal boosting obtained by low-pass filtering a point process.Darker shades indicate higher values in the contour plots, and the range is shared for the filtered processes.

Figure 3 :
Figure 3: The application of filtering to hickory and maple tress from the Lansing Woods data.The points, estimated spectral density function and low-pass filtered processes are shown from top to bottom respectively.In the middle plots, the shaded region delimited by the vertical line indicates the wavenumbers included in the filter.Darker shades in the contour plot indicate larger values, and the range is shared.

Figure 4 :
Figure 4: Some possible regions in wavenumber space which generalize a band.The transfer function is shown on the top row (where the function is either one, shown in gray, or zero, shown in white), and the corresponding impulse response on the bottom, with darker colors indicating larger values.To improve visualization, the impulse response functions are shown on a pseudo-log scale (with base 10).