NaturalNeighbours

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.

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.

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.

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.

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.

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]).

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.

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.

Sibson(d=0)

Interpolate using Sibson's coordinates with C(d) continuity at the data sites.

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).

Farin()

Interpolate using Farin's C(1) interpolant.

Laplace()

Interpolate using Laplace's coordinates.

Triangle(; allow_cache = true)

Interpolate using a piecewise linear interpolant over the underlying triangulation.

Nearest()

Interpolate by taking the function value at the nearest data site.

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.

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.

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.