Struct 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. The gradients are estimated with the estimate_gradients function. With flatness = 0.5 With flatness = 1.0

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 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 (should be >= 0.0) blends between an interpolation that adheres less to the gradient (small values) or adheres to it strongly (values larger than ~2.0) in the vicinity of any vertex. A value of 0.5 reproduces the original Sibson C1 interpolant. A value of 1.0 is the fastest and works well in many situations. A value of 0.0 reproduces the I1 interpolant introduced in Flötotto’s thesis chapter 4.4. If both the flatness and the gradients are 0.0 Sibsons C0 interpolant is recreated and Self::interpolate should be used.

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

Refer to NaturalNeighbor for more information and a visual example.

Example 1

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

struct PointWithData {
    position: Point2<f64>,
    height: f64,
    grad: [f64; 2],
}

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

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

let nn = triangulation.natural_neighbor();

// Interpolate point at coordinates (1.0, 2.0). This example uses a known gradients at each vertex.
// The gradient is the direction of steepest ascent for the function at that point
let query_point = Point2::new(1.0, 2.0);
let value: f64 = nn.interpolate_gradient(|v| v.data().height, |v| v.data().grad, 1.0, query_point).unwrap();

/// # Example 2

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).
let query_point = Point2::new(1.0, 2.0);
// Let's interpolate the value with estimated gradients. The gradient is the direction of steepest ascent for the function at that point
let value: f64 = nn.interpolate_gradient(|v| v.data().height, |v| nn.estimate_gradient(v, |h| h.data().height), 1.0, query_point).unwrap();

// could also have gone and calculated all grads upfront like this, notice the s in gradients
let grads = nn.estimate_gradients(|v| v.data().height);
let value: f64 = nn.interpolate_gradient(|v| v.data().height, &grads, 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”

It is also worth looking at Julia Flötotto’s thesis A Coordinate System associated to a Point Cloud issued from a Manifold: Definition, Properties and Applications chapter 4.1 and 4.4

pub fn estimate_gradients<I>( &self, i: I, ) -> impl Fn(VertexHandle<'_, V, DE, UE, F>) -> [<V as HasPosition>::Scalar; 2]

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

Estimates and returns the gradient for a every vertex in this triangulation.

This assumes that the triangulation models some kind of height field. The gradient is calculated from the weighted and normalized average of the normals of all triangles adjacent to the given vertex (weighted by triangle size).

NOTE: This is not Sibson’s gradient estimator

pub fn estimate_gradient<I>( &self, v: VertexHandle<'_, V, DE, UE, F>, i: I, ) -> [<V as HasPosition>::Scalar; 2]

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

Estimates and returns the gradient for a single vertex in this triangulation.

This assumes that the triangulation models some kind of height field. The gradient is calculated from the weighted and normalized average of the normals of all triangles adjacent to the given vertex (weighted by triangle size).

NOTE: This is not Sibson’s gradient estimator