NaturalNeighbours

Documentation for NaturalNeighbours.

This is a package for working with natural neighbours of planar point sets, enabling scattered data (or spatial) interpolation via natural neighbour interpolation, as well as derivative generation. We use DelaunayTriangulation.jl to define the spatial information. Much of the work in this package is based on this great thesis.

Please see the sections in the sidebar for some more discussion. The relevant docstrings for this package are shown below.

NaturalNeighbours.interpolateFunction
interpolate(tri::Triangulation, z; gradient=nothing, hessian=nothing, derivatives=false, kwargs...)
interpolate(points, z; gradient=nothing, hessian=nothing, derivatives=false, kwargs...)
interpolate(x::AbstractVector, y::AbstractVector, z; gradient=nothing, hessian=nothing, derivatives=false, kwargs...)

Construct an interpolant over the data z at the sites defined by the triangulation tri (or points, or (x, y)). See the Output section for a description of how to use the interpolant itp.

Missing vertices

When the underlying triangulation, tri, has points in get_points(tri) that are not vertices of the triangulation itself, the associated derivatives (relevant only if derivatives=true) at these points will be set to zero.

Keyword Arguments

  • gradient=nothing: The gradients at the corresponding data sites of z. Will be generated if isnothing(gradient) and derivatives==true.
  • hessian=nothing: The hessians at the corresponding data sites of z. Will be generated if isnothing(hessian) and derivatives==true.
  • derivatives=false: Whether to generate derivatives at the data sites of z. See also generate_derivatives.
  • kwargs...: Keyword arguments passed to generate_derivatives.

Output

The returned value is a NaturalNeighboursInterpolant struct. This struct is callable, with the following methods defined:

(itp::NaturalNeighboursInterpolant)(x, y, id::Integer=1; parallel=false, method=Sibson(), project = true, kwargs...)
(itp::NaturalNeighboursInterpolant)(vals::AbstractVector, x::AbstractVector, y::AbstractVector; parallel=true, method=Sibson(), project = true, kwargs...)
(itp::NaturalNeighboursInterpolant)(x::AbstractVector, y::AbstractVector; parallel=true, method=Sibson(), project = true, kwargs...)
  1. The first method is for scalars, with id referring to a thread id.
  2. This method is an in-place method for vectors, storing itp(x[i], y[i]) into vals[i].
  3. This method is similar to (2), but vals is constructed and returned.

In each method, method defines the method used for evaluating the interpolant, which is some AbstractInterpolator. For the first method, parallel is ignored, but for the latter two methods it defines whether to use multithreading or not for evaluating the interpolant at all the points. The kwargs... argument is passed into add_point! from DelaunayTriangulation.jl, e.g. you could pass some rng. Lastly, the project argument determines whether extrapolation is performed by projecting any exterior points onto the boundary of the convex hull of the data sites and performing two-point interpolation, or to simply replace any extrapolated values with Inf.

Performance

For the best performance when evaluating the interpolant at many points, either of the second or third methods are preferred over repeatedly calling the first.

Warning

Until we implement ghost point extrapolation, behaviour near the convex hull of your data sites may in some cases be undesirable, despite the extrapolation method we describe above, even for points that are inside the convex hull. If you want to control this behaviour so that you discard any points that are very close to the convex hull, see identify_exterior_points and the tol keyword argument.

source
NaturalNeighbours.differentiateFunction
differentiate(itp::NaturalNeighboursInterpolant, order)

Differentiate a given interpolant itp up to degree order (1 or 2). The returned object is a NaturalNeighboursDifferentiator struct, which is callable.

Missing vertices

When the underlying triangulation, tri, has points in get_points(tri) that are not vertices of the triangulation itself, the associated derivatives at these points will be set to zero.

For calling the resulting struct, we define the following methods:

(∂::NaturalNeighboursDifferentiator)(x, y, zᵢ, nc, id::Integer=1; parallel=false, method=default_diff_method(∂), kwargs...)
(∂::NaturalNeighboursDifferentiator)(x, y, id::Integer=1; parallel=false, method=default_diff_method(∂), interpolant_method=Sibson(), rng=Random.default_rng(), project = true, kwargs...)
(∂::NaturalNeighboursDifferentiator)(vals::AbstractVector, x::AbstractVector, y::AbstractVector; parallel=true, method=default_diff_method(∂), interpolant_method=Sibson(), kwargs...)
(∂::NaturalNeighboursDifferentiator{I, O})(x::AbstractVector, y::AbstractVector; parallel=true, method=default_diff_method(∂), interpolant_method=Sibson(), kwargs...) where {I, O}
  1. This method is useful if you already have an estimate for the function value, zᵢ, at the data site, (x, y), provided you also provide the NaturalCoordinates nc. id is the thread id.
  2. This method is for scalars, with id referring to a thread id.
  3. This method is an in-place method for vectors, storing ∂(x[i], y[i], 1) into vals[i].
  4. This method is similar to (3), but vals is constructed and returned.

The available keyword arguments are:

  • parallel=true: Whether to use multithreading. Ignored for the first two methods.
  • method=default_diff_method(∂): Default method for evaluating the interpolant. default_diff_method(∂) returns Direct(). The method must be a AbstractDifferentiator.
  • interpolant_method=Sibson(): The method used for evaluating the interpolant to estimate zᵢ for the latter three methods. See AbstractInterpolator for the available methods.
  • rng=Random.default_rng(): The random number generator used for estimating zᵢ for the latter three methods, or for constructing the natural coordinates.
  • project=false: Whether to project any extrapolated points onto the boundary of the convex hull of the data sites and perform two-point interpolation, or to simply replace any extrapolated values with Inf, when evaluating the interpolant in the latter three methods.
  • use_cubic_terms=true: If estimating second order derivatives, whether to use cubic terms. Only relevant for method == Direct().
  • alpha=0.1: If estimating second order derivatives, the weighting parameter used for estimating the second order derivatives. Only relevant for method == Iterative().
  • use_sibson_weight=true: Whether to weight the residuals in the associated least squares problems by the associated Sibson coordinates. Only relevant for method == Iterative() if order == 2.

The outputs are:

  • order == 1: The scalar methods return a Tuple of the form (∂x, ∂y), while the vector methods return a vector of Tuples of the form (∂x, ∂y).
  • order == 2: The scalar methods return a (∇, ℋ), where is a Tuple of the form (∂x, ∂y) and is a Tuple of the form (∂xx, ∂yy, ∂xy). The vector methods return a vector of (∇, ℋ)s.
Warning

Until we implement ghost point extrapolation, behaviour near the convex hull of your data sites may in some cases be undesirable, despite the extrapolation method we describe above, even for points that are inside the convex hull. If you want to control this behaviour so that you discard any points that are very close to the convex hull, see identify_exterior_points and the tol keyword argument.

source
NaturalNeighbours.generate_gradientsFunction
generate_gradients(
    tri,
    z,
    derivative_caches=[DerivativeCache(tri) for _ in 1:Base.Threads.nthreads()],
    neighbour_caches=[NaturalNeighboursCache(tri) for _ in 1:Base.Threads.nthreads()];
    parallel=true
)

Generate gradients at the data sites defined by the triangulation tri with associated function values tri.

Arguments

  • tri: A Triangulation object.
  • z: A vector of function values at the data sites.
  • derivative_caches=[DerivativeCache(tri) for _ in 1:Base.Threads.nthreads()]: A vector of DerivativeCache objects, one for each thread.
  • neighbour_caches=[NaturalNeighboursCache(tri) for _ in 1:Base.Threads.nthreads()]: A vector of NaturalNeighboursCache objects, one for each thread.

Keyword Arguments

  • parallel=true: Whether to use multithreading or not.

Output

  • : A vector of gradients at the data sites. Each element is a Tuple defining the gradient entries.
source
NaturalNeighbours.generate_derivativesFunction
generate_derivatives(
    tri,
    z,
    derivative_caches=[DerivativeCache(tri) for _ in 1:Base.Threads.nthreads()],
    neighbour_caches=[NaturalNeighboursCache(tri) for _ in 1:Base.Threads.nthreads()];
    parallel=true,
    method=Direct(),
    use_cubic_terms=true,
    alpha=0.1,
    initial_gradients=dwrap(method) == Direct() ? nothing : generate_gradients(tri, z, derivative_caches, neighbour_caches; method=dwrap(method), parallel, rng)
)

Generate derivatives at the data sites defined by the triangulation tri with associated function values tri.

Arguments

  • tri: A Triangulation object.
  • z: A vector of function values at the data sites.
  • derivative_caches=[DerivativeCache(tri) for _ in 1:Base.Threads.nthreads()]: A vector of DerivativeCache objects, one for each thread.
  • neighbour_caches=[NaturalNeighboursCache(tri) for _ in 1:Base.Threads.nthreads()]: A vector of NaturalNeighboursCache objects, one for each thread.

Keyword Arguments

  • parallel=true: Whether to use multithreading or not.
  • method=Direct(): The method used for generating the derivatives. See AbstractDifferentiator.
  • use_cubic_terms=true: Whether to use cubic terms for estimating the second order derivatives. Only relevant for method == Direct().
  • alpha=0.1: The weighting parameter used for estimating the second order derivatives. Only relevant for method == Iterative().
  • initial_gradients=dwrap(method) == Direct() ? nothing : generate_gradients(tri, z, derivative_caches, neighbour_caches; method=dwrap(method), parallel): The initial gradients used for estimating the second order derivatives. Only relevant for method == Iterative().

Output

  • : A vector of gradients at the data sites. Each element is a Tuple defining the gradient entries.
  • : A vector of Hessians at the data sites. Each element is a Tuple defining the Hessian entries in the form (H[1, 1], H[2, 2], H[1, 2]) (H[2, 1] is the same as H[2, 2]).
source
NaturalNeighbours.AbstractInterpolatorType
abstract type AbstractInterpolator{D}

Abstract type for defining the method used for evaluating an interpolant. D is, roughly, defined to be the smoothness at the data sites (currently only relevant for Sibson). The available subtypes are:

  • Sibson(d): Interpolate via the Sibson interpolant, with C(d) continuity at the data sites. Only defined for D ∈ (0, 1). If D == 1, gradients must be provided.
  • Triangle(d): Interpolate based on vertices of the triangle that the point resides in, with C(0) continuity at the data sites. D is ignored.
  • Nearest(d): Interpolate by returning the function value at the nearest data site. D doesn't mean much here (it could be D = ∞), and so it is ignored and replaced with 0.
  • Laplace(d): Interpolate via the Laplace interpolant, with C(0) continuity at the data sites. D is ignored.
  • Farin(d): Interpolate using the Farin interpolant, with C(1) continuity at the data sites. d is ignored.
  • Hiyoshi(d): Interpolate using the Hiyoshi interpolant, with C(d) continuity at the data sites. Currently, only defined for d == 2.

Our implementation of Sibson(0)'s coordinates follows this article with some simple modifications.

source
NaturalNeighbours.AbstractDifferentiatorType
abstract type AbstractDifferentiator end

Abstract type for defining the method used for differentiating an interpolant or generating derivatives at data sites.

  • Direct(): Generate derivatives directly with one least squares problem.
  • Iterative(): Generate derivatives iteratively: Gradients are estimated first, and then both gradients and Hessians are estimated with the initial gradients used to refine the results.
source
NaturalNeighbours.HiyoshiType
Hiyoshi(d)

Interpolate using Hiyoshi's C(d) interpolant. Hiyoshi's interpolant C(0) is not yet implemented, but we do not make any conversions to C(2) like in e.g. Farin(), e.g. Farin() gets converted to Farin(1) but, to support possible later versions, Hiyoshi() does not get converted to Hiyoshi(2).

source
NaturalNeighbours.TriangleType
Triangle(; allow_cache = true)

Interpolate using a piecewise linear interpolant over the underlying triangulation.

Cached coordinates with `allow_cache=true`

The Triangle() interpolator is special as it will cache the coordinates used for each triangle. In particular, when an interpolator is evaluated with the Triangle() method, the object returned from Triangle() will store all the coordinates. For this reason, if you want to reuse Triangle() for different evaluations of the interpolant, you should be sure to reuse the same instance rather than reinstantiating it every single time. If you do not want this behaviour, set allow_cache = false.

If you only ever call the scalar-argument version of the interpolant, no caching will be done even with allow_cache = true.

source
NaturalNeighbours.IterativeType
Iterative()

Generate derivatives iteratively: Gradients are estimated first, and then both gradients and Hessians are estimated with the initial gradients used to refine the results.

source
NaturalNeighbours.identify_exterior_pointsFunction
identify_exterior_points(x, y, points, boundary_nodes; tol = 0.0)

Given a polygon described by (points, boundary_nodes), matching the specification of polygons in DelaunayTriangulation.jl, returns a vector of indices of the points defined by (x, y) that are outside of the polygon.

Use tol to specify a tolerance for the distance to the polygon.

source
identify_exterior_points(x, y, itp::NaturalNeighboursInterpolant; tol = 0.0)

Returns the indices of the points defined by the vectors (x, y) that are outside of the underlying triangulation to the interpolant itp.

Use tol to specify a tolerance for the distance to the triangulation.

source