Struct spade::NaturalNeighbor

pub struct NaturalNeighbor<'a, T>
where
    T: Triangulation,
    T::Vertex: HasPosition,
{ /* private fields */ }

Implements methods for natural neighbor interpolation.

Natural neighbor interpolation is a spatial interpolation method. For a given set of 2D input points with an associated value (e.g. a height field of some terrain or temperature measurements at different locations in a country), natural neighbor interpolation allows to smoothly interpolate the associated value for every location within the convex hull of the input points.

Spade currently assists with 4 interpolation strategies:

• Nearest neighbor interpolation: Fastest. Exhibits too poor quality for many tasks. Not continuous along the edges of the voronoi diagram. Use DelaunayTriangulation::nearest_neighbor.
• Barycentric interpolation: Fast. Not smooth on the edges of the Delaunay triangulation. See Barycentric.
• Natural neighbor interpolation: Slower. Smooth everywhere except the input points. The input points have no derivative. See NaturalNeighbor::interpolate
• Natural neighbor interpolation with gradients: Slowest. Smooth everywhere, even at the input points. See NaturalNeighbor::interpolate_gradient.

Performance comparison

In general, the speed of interpolating random points in a triangulation will be determined by the speed of the lookup operation after a certain triangulation size is reached. Thus, if random sampling is required, using a crate::HierarchyHintGenerator may be beneficial.

If subsequent queries are close to each other, crate::LastUsedVertexHintGenerator should also work fine. In this case, interpolating a single value should result in the following (relative) run times:

nearest neighbor	barycentric	natural neighbor (no gradients)	natural neighbor (gradients)
23ns	103ns	307ns	422ns

Usage

This type is created by calling DelaunayTriangulation::natural_neighbor. It contains a few internal buffers that are used to prevent recurring allocations. For best performance it should be created only once per thread and then used in all interpolation activities (see example).

Example

use spade::{Point2, HasPosition, DelaunayTriangulation, InsertionError, Triangulation as _};

struct PointWithHeight {
    position: Point2<f64>,
    height: f64,
}

impl HasPosition for PointWithHeight {
    type Scalar = f64;

    fn position(&self) -> Point2<f64> {
        self.position
    }
}

fn main() -> Result<(), InsertionError> {
    let mut t = DelaunayTriangulation::<PointWithHeight>::new();
    t.insert(PointWithHeight { position: Point2::new(-1.0, -1.0), height: 42.0})?;
    t.insert(PointWithHeight { position: Point2::new(-1.0, 1.0), height: 13.37})?;
    t.insert(PointWithHeight { position: Point2::new(1.0, -1.0), height: 27.18})?;
    t.insert(PointWithHeight { position: Point2::new(1.0, 1.0), height: 31.41})?;

    // Set of query points (would be many more realistically):
    let query_points = [Point2::new(0.0, 0.1), Point2::new(0.5, 0.1)];

    // Good: Re-use interpolation object!
    let nn = t.natural_neighbor();
    for p in &query_points {
        println!("Value at {:?}: {:?}", p, nn.interpolate(|v| v.data().height, *p));
    }

    // Bad (slower): Don't re-use internal buffers
    for p in &query_points {
        println!("Value at {:?}: {:?}", p, t.natural_neighbor().interpolate(|v| v.data().height, *p));
    }
    Ok(())
}

Visual comparison of interpolation algorithms

Note: All of these images are generated by the “interpolation” example

Nearest neighbor interpolation exhibits discontinuities along the voronoi edges:

Barycentric interpolation, by contrast, is continuous everywhere but has no derivative on the edges of the Delaunay triangulation. These show up as sharp corners in the color gradients on the drawn edges:

By contrast, natural neighbor interpolation is smooth on the edges - the previously sharp angles are now rounded off. However, the vertices themselves are still not continuous and will form sharp “peaks” in the resulting surface.

With a gradient, the sharp peaks are gone - the surface will smoothly approximate a linear function as defined by the gradient in the vicinity of each vertex. In the image below, a gradient of (0.0, 0.0) is used which leads to a small “disc” around each vertex with values close to the vertex value.

Implementations

impl<'a, V, DE, UE, F, L> NaturalNeighbor<'a, DelaunayTriangulation<V, DE, UE, F, L>>

where V: HasPosition, DE: Default, UE: Default, F: Default, L: HintGenerator<<V as HasPosition>::Scalar>, <V as HasPosition>::Scalar: Float,

pub fn get_weights( &self, position: Point2<<V as HasPosition>::Scalar>, result: &mut Vec<(FixedVertexHandle, <V as HasPosition>::Scalar)> )

Calculates the natural neighbors and their weights (sibson coordinates) of a given query position.

The neighbors are returned in clockwise order. The weights will add up to 1.0. The neighbors are stored in the result parameter to prevent unnecessary allocations. result will be cleared initially.

The number of returned natural neighbors depends on the given query position:

• result will be empty if the query position lies outside of the triangulation’s convex hull
• result will contain exactly one vertex if the query position is equal to that vertex position.
• result will contain exactly two entries if the query position lies exactly on an edge of the convex hull.
• result will contain at least three (vertex, weight) tuples if the query point lies on an inner face or an inner edge.

Example: The natural neighbors (red vertices) of the query point (blue dot) with their weights. The elements will be returned in clockwise order as indicated by the indices drawn within the red circles.

pub fn interpolate<I>( &self, i: I, position: Point2<<V as HasPosition>::Scalar> ) -> Option<<V as HasPosition>::Scalar>

where I: Fn(VertexHandle<'_, V, DE, UE, F>) -> <V as HasPosition>::Scalar,

Interpolates a value at a given position.

Returns None for any point outside the triangulations convex hull. The value to interpolate is given by the i parameter. The resulting interpolation will be smooth everywhere except at the input vertices.

Refer to NaturalNeighbor for an example on how to use this function.

pub fn interpolate_gradient<I, G>( &self, i: I, g: G, flatness: <V as HasPosition>::Scalar, position: Point2<<V as HasPosition>::Scalar> ) -> Option<<V as HasPosition>::Scalar>

where I: Fn(VertexHandle<'_, V, DE, UE, F>) -> <V as HasPosition>::Scalar, G: Fn(VertexHandle<'_, V, DE, UE, F>) -> [<V as HasPosition>::Scalar; 2],

Interpolates a value at a given position.

In contrast to Self::interpolate, this method has a well defined derivative at each vertex and will approximate a linear function in the proximity of any vertex.

The value to interpolate is given by the i parameter. The gradient that defines the derivative at each input vertex is given by the g parameter.

The flatness parameter blends between an interpolation that ignores the given gradients (value 0.0) or adheres to it strongly (values larger than ~2.0) in the vicinity of any vertex. When in doubt, using a value of 1.0 should result in a good interpolation and is also the fastest.

Returns None for any point outside of the triangulation’s convex hull.

Refer to NaturalNeighbor for more information and a visual example.

Example

use spade::{DelaunayTriangulation, HasPosition, Point2};

struct PointWithHeight {
    position: Point2<f64>,
    height: f64,
}

impl HasPosition for PointWithHeight {
    type Scalar = f64;
    fn position(&self) -> Point2<f64> { self.position }
}

let mut triangulation: DelaunayTriangulation<PointWithHeight> = Default::default();
// Insert some points into the triangulation
triangulation.insert(PointWithHeight { position: Point2::new(10.0, 10.0), height: 0.0 });
triangulation.insert(PointWithHeight { position: Point2::new(10.0, -10.0), height: 0.0 });
triangulation.insert(PointWithHeight { position: Point2::new(-10.0, 10.0), height: 0.0 });
triangulation.insert(PointWithHeight { position: Point2::new(-10.0, -10.0), height: 0.0 });

let nn = triangulation.natural_neighbor();

// Interpolate point at coordinates (1.0, 2.0). This example uses a fixed gradient of (0.0, 0.0) which
// means that the interpolation will have normal vector parallel to the z-axis at each input point.
// Realistically, the gradient might be stored as an additional property of `PointWithHeight`.
let query_point = Point2::new(1.0, 2.0);
let value: f64 = nn.interpolate_gradient(|v| v.data().height, |_| [0.0, 0.0], 1.0, query_point).unwrap();

References

This method uses the C1 extension proposed by Sibson in “A brief description of natural neighbor interpolation, R. Sibson, 1981”