This command computes a diffusion kernel estimate
of point process intensity from the observed point pattern `x`

.

The function `densityHeat`

is generic,
with methods for point patterns in two dimensions
(class `"ppp"`

) and point patterns on a linear network
(class `"lpp"`

). The function `densityHeat.ppp`

described
here is the method for class `"ppp"`

. Given a two-dimensional
point pattern `x`

, it computes a diffusion kernel estimate
of the intensity of the point process which generated `x`

.

Diffusion kernel estimates were developed
by Botev et al (2010), Barry and McIntyre (2011) and
Baddeley et al (2021).

Barry and McIntyre (2011) proposed an estimator for point process
intensity based on a random walk on the pixel grid inside the
observation window. Baddeley et al (2021) showed that the
Barry-McIntyre method is a special case of the *diffusion estimator*
proposed by Botev et al (2010).

The original Barry-McIntyre algorithm assumes a symmetric random walk
(i.e. each possible transition has the same probability \(p\))
and requires a square pixel grid (i.e. equal
spacing in the \(x\) and \(y\) directions). Their original
algorithm is used if `symmetric=TRUE`

. Use the `...`

arguments to ensure a square grid: for example, the argument
`eps`

specifies a square grid with spacing `eps`

units.

The more general algorithm used here (Baddeley et al, 2021)
does not require a square grid of pixels.
If the pixel grid is not square, and if `symmetric=FALSE`

(the default), then the random walk is not symmetric,
in the sense that the probabilities of different jumps will be
different, in order to ensure that the smoothing is isotropic.

This implementation also includes two generalizations to
the case of adaptive smoothing (Baddeley et al, 2021).

In the first version of adaptive smoothing, the bandwidth is
spatially-varying.
The argument `sigma`

should be a pixel image (class `"im"`

)
or a `function(x,y)`

specifying the bandwidth at each spatial
location. The smoothing is performed by solving the
heat equation with spatially-varying parameters.

In the second version of adaptive smoothing, each data point in
`x`

is smoothed using a separate bandwidth.
The argument `sigmaX`

should be a numeric vector
specifying the bandwidth for each point of `x`

.
The smoothing is performed using the lagged arrival algorithm.
The argument `sigma`

can be omitted.

If `extrapolate=FALSE`

(the default), calculations are performed
using the Euler scheme for the heat equation.
If `extrapolate=TRUE`

, the accuracy of the result will be
improved by applying Richardson extrapolation (Baddeley et al, 2021, Section
4). After computing the intensity estimate using the Euler scheme
on the desired pixel grid, another estimate is computed using the same
method on another pixel grid, and the two estimates are combined by
Richardson extrapolation to obtain a more accurate result.
The second grid is coarser than the original grid if
`coarsen=TRUE`

(the default), and finer than the original grid
if `coarsen=FALSE`

. Setting `extrapolate=TRUE`

increases
computation time by 35% if `coarsen=TRUE`

and by 400% if
`coarsen=FALSE`

.