Hi all,
I’ve just submitted my package NaturalNeighbours.jl for registration. This package is useful for performing natural neighbour interpolation over (scattered) planar data sets, as well as for derivative generation. The docs hopefully define everything in sufficient detail. I summarise some key parts of the package below too.
The idea behind natural neighbour interpolation is simple. Given some data (\boldsymbol x_i, z_i), i=1,\ldots,m, where \boldsymbol x_i \in \boldsymbol X \subseteq \mathbb R^2, and z_i \in Z \subseteq \mathbb R, suppose we want to interpolate the data at some \boldsymbol x_0. The idea is to use spatial information to decide which points to use for estimating the function at \boldsymbol x_0, rather than having to think hard about features like the number of nearest neighbours to use or a cut-off radius common in other spatial interpolation methods. Spatial information is represented via the Voronoi tessellation of the data sites, constructed with my other package DelaunayTriangulation.jl, allowing for the definition of the natural neighbours to \boldsymbol x_0, denoted N_0, to be the set of points in \boldsymbol X whose Voronoi cell contains \boldsymbol x_0. Then, using some definition of local coordinates \boldsymbol \lambda depending on \boldsymbol x_0 and N_0, the interpolant we define is
Thanks to certain properties of the Voronoi tessellation, this interpolant is continuous in the query point \boldsymbol x_0 and smooth (only having derivative discontinuities at the data sites). Moreover, since this is a local interpolant, it is extremely fast to evaluate – I can evaluate the interpolant at hundreds of thousands of data points in less than half a second on my machine, independent of the size of the original data set. More mathematical detail is available here, where I also define the available coordinates: Sibson’s, Sibson-1’s (differentiable at the data sites!), Triangle, Nearest, and Laplace. A detailed example demonstrating how we can use interpolation is given here (and another example is at the end of this post).
Methods for derivative generation are also available – generation is used rather than the term estimation since the goal of these methods is to reconstruct derivatives that lead to a satisfactory interpolant, rather than necessarily reconstructing the exact values of the underlying function (see this paper for more discussion on this point if you are interested). By defining certain weighted least squares problems, gradients and Hessians can be accurately computed at the provided data sites, and the interpolants we provide can also be differentiated. The mathematics behind this is given here, and a detailed example of differentiation is given here.
A quick demonstration of the package is below.
using NaturalNeighbours
using CairoMakie
using StableRNGs
## The data
rng = StableRNG(123)
f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)
x = vec([(i - 1) / 9 for i in (1, 3, 4, 5, 8, 9, 10), j in (1, 2, 3, 5, 6, 7, 9, 10)])
y = vec([(j - 1) / 9 for i in (1, 3, 4, 5, 8, 9, 10), j in (1, 2, 3, 5, 6, 7, 9, 10)])
z = f.(x, y)
## The interpolant and grid
itp = interpolate(x, y, z; derivatives=true)
xg = LinRange(0, 1, 100)
yg = LinRange(0, 1, 100)
_x = vec([x for x in xg, _ in yg])
_y = vec([y for _ in xg, y in yg])
exact = f.(_x, _y)
## Evaluate some interpolants
sibson_vals = itp(_x, _y; method=Sibson())
triangle_vals = itp(_x, _y; method=Triangle())
laplace_vals = itp(_x, _y; method=Laplace())
sibson_1_vals = itp(_x, _y; method=Sibson(1))
nearest = itp(_x, _y; method=Nearest())
## Plot
function plot_2d(fig, i, j, title, vals, xg, yg, x, y, show_scatter=true)
ax = Axis(fig[i, j], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
contourf!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
show_scatter && scatter!(ax, x, y, color=:red, markersize=14)
end
function plot_3d(fig, i, j, title, vals, xg, yg)
ax = Axis3(fig[i, j], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
surface!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
end
all_vals = (sibson_vals, triangle_vals, laplace_vals, sibson_1_vals, nearest, exact)
titles = ("(a): Sibson", "(b): Triangle", "(c): Laplace", "(d): Sibson-1", "(e): Nearest", "(f): Exact")
fig = Figure(fontsize=36)
for (i, (vals, title)) in enumerate(zip(all_vals, titles))
plot_2d(fig, 1, i, title, vals, xg, yg, x, y, !(vals === exact))
plot_3d(fig, 2, i, " ", vals, xg, yg)
end
resize_to_layout!(fig)
fig
# could keep going and differentiating, etc...
# ∂ = differentiate(itp, 2) -- see the docs.
I also give an example of using interpolation for studying Switzerland’s elevation data here. A figure from this example is shown below.
I don’t see any errors. Maybe it’s a browser issue? Unless I’m just missing them
