GitHub

Interpolation

For our first application, we consider applying Voronoi tessellations to interpolation. The discussion here is based on the implementation of what is known as natural neighbour interpolation in NaturalNeighbours.jl. NaturalNeighbours.jl considers a diverse set of interpolants for this purpose, but here we will focus only on a single interpolant (called the Sibson() interpolant in NaturalNeighbours.jl).

The problem we are considering is as follows: Given some data $\mathcal D = \{(\vb x_i, z_i)\}_{i=1}^m \subseteq \mathbb R^2 \times \mathbb R$, we want to construct a smooth interpolant $f \colon \mathcal C\mathcal H(X) \to \mathbb R$ of the data, where $X = \{\vb x_1, \ldots, \vb x_m\}$, such that $f(\vb x_i) = z_i$ for all $i = 1, \ldots, m$.

One idea to interpolate this data would be to use a piecewise linear interpolant, obtained by defining, inside each triangle $T \in \mathcal D\mathcal T(X)$, a piecewise linear interpolant using the data points at each of the three vertices. One problem with this is that $\mathcal D\mathcal T(X)$ does not depend continuously on $X$, meaning that small perturbations of a data point can lead to topological changes in the triangulation. One way around this is to instead use the Voronoi tessellation to guide the tessellation, since $\mathcal V(X)$ does depend continuously on $X$. We give an example of this below, where we show how $\mathcal D\mathcal T(X)$ may change significant after a small perturbation while $\mathcal V(X)$ does not at the same time.

Example block output

This observation motivates the use of the Voronoi tessellation to guide the interpolation.

Natural Neighbours

Let us start by defining natural neighbours. Consider two Voronoi polygons $\mathcal V_i$ and $\mathcal V_j$, and let $\mathcal F_{ij} = \mathcal V_i \cap \mathcal V_j$. If $\mathcal F_{ij} \neq \emptyset$, we say that $\vb x_i$ and $\vb x_j$ are natural neighbours in $X$. For a point $\vb \in X$, we denote its set of natural neighbours by $N(\vb x) \subseteq X$, and the corresponding indices by $N_i = \{j : \vb x_j \in N(\vb x_i)\}$.

We use natural neighbours to guide our representation of points. In particular, we use natural neighbour coordinates, which are defined as follows: For a point $\vb x \in X$, a vector $\boldsymbol\lambda$ are called a set of natural neighbour coordinates if $\lambda_{i} \geq 0$ for each $i$, $\boldsymbol\lambda$ is continuous with respect to $\vb x$, and $\lambda_i > 0 \iff \vb x_i \in N(\vb x)$.

Sibsonian Interpolation

Here, for $\boldsymbol\lambda$ we will discuss Sibsonian coordinates. To define these coordinates, take $\mathcal V(X)$ and consider what happens when a new point $\vb x_0$ is added into it. This new point creates a new tessellation $\mathcal V(X \cup \{\vb x_0\})$, where all the tiles that change in the new tessellation compared to $\mathcal V(\vb X)$ are those in $N(\vb x_0)$. We let $A(\vb x_0)$ be the area of the new tile created by $\vb x_0$, and let $A(\vb x_i)$ be the area of the intersection between the original tile for $\vb x_i$ and the new tile from $\vb x_0$. The Sibson coordinates are then

\[\lambda_i(\vb x_0) = \frac{A(\vb x_i)}{A(\vb x_0)}.\]

Using this definition of Sibsonian coordinates, the Sibson interpolant is defined by

\[f(\vb x_0) = \sum_{i \in N_0} \lambda_i(\vb x_0)z_i.\]

This interpolant is $C^1$ continuous in $\mathcal C\mathcal H(X) \setminus X$, with derivative discontinuities at the data sites.

Implementation

Let us now implement the Sibson interpolant. We note that, in this implementation, we ignore some edge cases; these are handled properly in the implementation in NaturalNeighbours.jl. Moreover, our implementation will not be the most efficient, but will be enough for the purposes of this demonstration. We assume that all points are contained in the interior of $\mathcal C\mathcal H(X)$.

Bowyer-Watson envelope

The first issue to deal with in our implementation is the computation of the Bowyer-Watson envelope. For a point $\vb x$, its Bowyer-Watson envelope is the boundary of the set of all triangles in $\mathcal D\mathcal T(X)$ whose circumcircles contain $\vb x$. This envelope tells us the region in which any changes to the triangulation and to the tessellation can occur. A reasonably straight forward way to implement this is to simply add $\vb x$ into $\mathcal D\mathcal T(X)$ and take all the triangles containing $\vb x$ as a vertex in $\mathcal D\mathcal T(X \cup \{\vb x\})$. We then remove $\vb x$ from $\mathcal D\mathcal T(X \cup \{\vb x\})$. We implement this below.[1]

using DelaunayTriangulation
function compute_envelope(tri::Triangulation, point)
    r = DelaunayTriangulation.num_points(tri)
    x, y = getxy(point)
    add_point!(tri, x, y)
    envelope_vertices = DelaunayTriangulation.get_surrounding_polygon(tri, r + 1)
    push!(envelope_vertices, envelope_vertices[begin])
    envelope_points = [get_point(tri, i) for i in envelope_vertices]
    delete_point!(tri, r + 1)
    DelaunayTriangulation.pop_point!(tri)
    return envelope_vertices, envelope_points
end
compute_envelope (generic function with 1 method)

Let's now check that this function works.

using CairoMakie
using StableRNGs
using ElasticArrays
rng = StableRNG(999)
points = ElasticMatrix(randn(rng, 2, 50)) # so that the points are mutable
tri = triangulate(points; rng)
envelope_vertices, envelope_points = compute_envelope(tri, (0.5, 0.5))

fig = Figure(fontsize = 24)
ax = Axis(fig[1, 1], width = 400, height = 400)
triplot!(ax, tri, show_points = true)
ax2 = Axis(fig[1, 2], width = 400, height = 400)
add_point!(tri, 0.5, 0.5)
triplot!(ax2, tri, show_points = true)
poly!(ax2, envelope_points, color = (:red, 0.2))
resize_to_layout!(fig)
fig
Example block output

As we can see, the red region we have computed from our envelope is indeed the envelope we need.

Computing the Sibsonian coordinates

When we are interpolating at a point $\vb x$, remember that we need to know the area of the Voronoi polygon that would be produced when $\vb x$ is inserted into $\mathcal V(X)$. To compute this area, we need to know how we can compute it using the Bowyer-Watson envelope. Remember that the Voronoi polygons are obtained by drawing lines between the circumcenters of neighbouring triangles. This can be done using the Bowyer-Watson envelope: The edges of the envelope, together with the new point, define the triangles that would be produced if it were to be added into the triangulation, and so we can join the circumcenters of these triangles to compute the new Voronoi polygon. We can then use this polygon, together with the original tessellation, to compute the Sibsonian coordinates.

Let us first discuss the area of the Voronoi polygons in $N(\vb x)$ ($\vb x$'s natural neighbours, i.e. the vertices of the envelope) before $\vb x$ is inserted. We only need to compute the part of the area that is contained within the envelope, since everything outside of that envelope is unchanged. If we included the entire area, then the area that we subtract off for the intersection we compute later would just cancel it out anyway. Let's zoom in on the envelope and consider a specific example of how we can do this computation.

fig = Figure(fontsize = 24)
ax = Axis(fig[1, 1], width = 400, height = 400)
triplot!(ax, tri, show_points = true)
lines!(ax, envelope_points, color = :red)
j = 7 # example vertex
v = envelope_vertices[j]
scatter!(ax, [get_point(tri, v)], color = :blue)
first_neighbour = envelope_vertices[j - 1]
next_triangle = get_adjacent(tri, first_neighbour, v)
next_triangle_2 = get_adjacent(tri, v, envelope_vertices[j + 1])
last_neighbour = envelope_vertices[j + 1]
polygon_points = [
    get_point(tri, v),
    (get_point(tri, v) .+ get_point(tri, last_neighbour)) ./ 2,
    DelaunayTriangulation.triangle_circumcenter(tri, (v, envelope_vertices[j + 1], next_triangle_2)),
    DelaunayTriangulation.triangle_circumcenter(tri, (first_neighbour, v, next_triangle)),
    (get_point(tri, v) .+ get_point(tri, first_neighbour)) ./ 2,
]
poly!(ax, polygon_points, color = (:blue, 0.5), strokecolor = :blue, strokewidth = 2)
xlims!(ax, -0.5, 1.4)
ylims!(ax, -0.15, 1.4)
resize_to_layout!(fig)
fig
Example block output

The relevant polygon is shown above in blue, associated with the generator shown by the blue point. We need to compute the area of this polygon. This is simple using the shoelace formula. Our implementation of this is given below.

function polygon_area(points) # this is the first formula in the "Other formulae" section of the above Wikipedia article
    n = DelaunayTriangulation.num_points(points)
    p, q, r, s = get_point(points, 1, 2, n, n - 1)
    px, py = getxy(p)
    qx, qy = getxy(q)
    rx, ry = getxy(r)
    sx, sy = getxy(s)
    area = px * (qy - ry) + rx * (py - sy)
    for i in 2:(n - 1)
        p, q, r = get_point(points, i, i + 1, i - 1)
        px, py = getxy(p)
        qx, qy = getxy(q)
        rx, ry = getxy(r)
        area += px * (qy - ry)
    end
    return area / 2
end
function pre_insertion_area(tri::Triangulation, i, envelope_vertices) # area from the envelope[i]th generator
    poly_points = NTuple{2, Float64}[]
    u = envelope_vertices[i]
    prev_index = i == 1 ? length(envelope_vertices) - 1 : i - 1
    next_index = i == length(envelope_vertices) ? 1 : i + 1
    first_neighbour = envelope_vertices[prev_index]
    last_neighbour = envelope_vertices[next_index]
    v = last_neighbour
    (ux, uy), (vx, vy) = get_point(tri, u, v)
    mx1, my1 = (ux + vx) / 2, (uy + vy) / 2
    push!(poly_points, (mx1, my1))
    while v ≠ first_neighbour
        w = get_adjacent(tri, u, v)
        cx, cy = DelaunayTriangulation.triangle_circumcenter(tri, (u, v, w))
        push!(poly_points, (cx, cy))
        v = w
    end
    vx, vy = get_point(tri, v)
    mx, my = (ux + vx) / 2, (uy + vy) / 2
    push!(poly_points, (mx, my), (mx1, my1))
    return polygon_area(poly_points)
end
pre_insertion_area (generic function with 1 method)

The details for the post-insertion area are similar, but now the triangles that we take the circumcenters of are those where the edges instead join with the inserted vertex. The function we use is below.

function post_insertion_area(tri::Triangulation, i, envelope_vertices, point)
    u = envelope_vertices[i]
    prev_index = i == 1 ? length(envelope_vertices) - 1 : i - 1
    next_index = i == length(envelope_vertices) ? 1 : i + 1
    first_neighbour = envelope_vertices[prev_index]
    last_neighbour = envelope_vertices[next_index]
    p, q, r = get_point(tri, u, first_neighbour, last_neighbour)
    px, py = getxy(p)
    qx, qy = getxy(q)
    rx, ry = getxy(r)
    mpq = (px + qx) / 2, (py + qy) / 2
    mpr = (px + rx) / 2, (py + ry) / 2
    g1 = DelaunayTriangulation.triangle_circumcenter(p, r, point)
    !all(isfinite, g1) && return NaN # point is one of p and r, i.e. we are interpolating at a data site
    g2 = DelaunayTriangulation.triangle_circumcenter(q, p, point)
    !all(isfinite, g2) && return NaN
    points = (mpq, mpr, g1, g2, mpq)
    return polygon_area(points)
end
post_insertion_area (generic function with 1 method)

Now that we can compute the pre- and post-insertion areas, we can start computing the Sibsonian coordinates.

function compute_sibson_coordinates(tri::Triangulation, envelope_vertices, point)
    coordinates = zeros(length(envelope_vertices) - 1)
    w = 0.0
    for i in firstindex(envelope_vertices):(lastindex(envelope_vertices) - 1)
        pre = max(0.0, pre_insertion_area(tri, i, envelope_vertices))
        post = max(0.0, post_insertion_area(tri, i, envelope_vertices, point))
        if isnan(post) # need to return the the vector λ = [1] since we are exactly at a data site
            return [1.0]
        end
        coordinates[i] = max(pre - post, 0.0) # take care of any precision issues
        w += coordinates[i]
    end
    coordinates ./= w
    return coordinates
end
compute_sibson_coordinates (generic function with 1 method)

This function gives our $\boldsymbol\lambda$ vector. Notice that, in the computation of these coordinates, we need needed to have $\mathcal V(X)$ directly or make use of the data $z_i$.

Evaluating the Sibsonian interpolant

Now we can evaluate our Sibson interpolant. The following function does this for us.

function evaluate_sibson_interpolant(tri::Triangulation, z, point)
    envelope_vertices, _ = compute_envelope(tri, point)
    λ = compute_sibson_coordinates(tri, envelope_vertices, point)
    if length(λ) == 1
        for i in each_solid_vertex(tri)
            get_point(tri, i) == point && return z[i]
        end
    else
        itp = 0.0
        for (λ, k) in zip(λ, envelope_vertices)
            itp += λ * z[k]
        end
        return itp
    end
end
evaluate_sibson_interpolant (generic function with 1 method)

Let's now use this function to interpolate some data.

f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)
trit = triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 30, 30)
zz = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(trit)]
xx = LinRange(0.001, 0.999, 20) # handling points on the boundary requires more care than we have discussed here
yy = LinRange(0.001, 0.999, 20)
fig = Figure(fontsize = 24)
ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"y", title = "True function", titlealign = :left, width = 400, height = 400)
contourf!(ax, xx, yy, f.(xx, yy'))
ax2 = Axis(fig[1, 2], xlabel = L"x", ylabel = L"y", title = "Interpolant", titlealign = :left, width = 400, height = 400)
zi = [evaluate_sibson_interpolant(trit, zz, (xᵢ, yᵢ)) for xᵢ in xx, yᵢ in yy]
contourf!(ax2, xx, yy, zi)
resize_to_layout!(fig)
fig
Example block output

Works perfectly!

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using DelaunayTriangulation
function compute_envelope(tri::Triangulation, point)
    r = DelaunayTriangulation.num_points(tri)
    x, y = getxy(point)
    add_point!(tri, x, y)
    envelope_vertices = DelaunayTriangulation.get_surrounding_polygon(tri, r + 1)
    push!(envelope_vertices, envelope_vertices[begin])
    envelope_points = [get_point(tri, i) for i in envelope_vertices]
    delete_point!(tri, r + 1)
    DelaunayTriangulation.pop_point!(tri)
    return envelope_vertices, envelope_points
end

using CairoMakie
using StableRNGs
using ElasticArrays
rng = StableRNG(999)
points = ElasticMatrix(randn(rng, 2, 50)) # so that the points are mutable
tri = triangulate(points; rng)
envelope_vertices, envelope_points = compute_envelope(tri, (0.5, 0.5))

fig = Figure(fontsize = 24)
ax = Axis(fig[1, 1], width = 400, height = 400)
triplot!(ax, tri, show_points = true)
ax2 = Axis(fig[1, 2], width = 400, height = 400)
add_point!(tri, 0.5, 0.5)
triplot!(ax2, tri, show_points = true)
poly!(ax2, envelope_points, color = (:red, 0.2))
resize_to_layout!(fig)
fig

fig = Figure(fontsize = 24)
ax = Axis(fig[1, 1], width = 400, height = 400)
triplot!(ax, tri, show_points = true)
lines!(ax, envelope_points, color = :red)
j = 7 # example vertex
v = envelope_vertices[j]
scatter!(ax, [get_point(tri, v)], color = :blue)
first_neighbour = envelope_vertices[j - 1]
next_triangle = get_adjacent(tri, first_neighbour, v)
next_triangle_2 = get_adjacent(tri, v, envelope_vertices[j + 1])
last_neighbour = envelope_vertices[j + 1]
polygon_points = [
    get_point(tri, v),
    (get_point(tri, v) .+ get_point(tri, last_neighbour)) ./ 2,
    DelaunayTriangulation.triangle_circumcenter(tri, (v, envelope_vertices[j + 1], next_triangle_2)),
    DelaunayTriangulation.triangle_circumcenter(tri, (first_neighbour, v, next_triangle)),
    (get_point(tri, v) .+ get_point(tri, first_neighbour)) ./ 2,
]
poly!(ax, polygon_points, color = (:blue, 0.5), strokecolor = :blue, strokewidth = 2)
xlims!(ax, -0.5, 1.4)
ylims!(ax, -0.15, 1.4)
resize_to_layout!(fig)
fig

function polygon_area(points) # this is the first formula in the "Other formulae" section of the above Wikipedia article
    n = DelaunayTriangulation.num_points(points)
    p, q, r, s = get_point(points, 1, 2, n, n - 1)
    px, py = getxy(p)
    qx, qy = getxy(q)
    rx, ry = getxy(r)
    sx, sy = getxy(s)
    area = px * (qy - ry) + rx * (py - sy)
    for i in 2:(n - 1)
        p, q, r = get_point(points, i, i + 1, i - 1)
        px, py = getxy(p)
        qx, qy = getxy(q)
        rx, ry = getxy(r)
        area += px * (qy - ry)
    end
    return area / 2
end
function pre_insertion_area(tri::Triangulation, i, envelope_vertices) # area from the envelope[i]th generator
    poly_points = NTuple{2, Float64}[]
    u = envelope_vertices[i]
    prev_index = i == 1 ? length(envelope_vertices) - 1 : i - 1
    next_index = i == length(envelope_vertices) ? 1 : i + 1
    first_neighbour = envelope_vertices[prev_index]
    last_neighbour = envelope_vertices[next_index]
    v = last_neighbour
    (ux, uy), (vx, vy) = get_point(tri, u, v)
    mx1, my1 = (ux + vx) / 2, (uy + vy) / 2
    push!(poly_points, (mx1, my1))
    while v ≠ first_neighbour
        w = get_adjacent(tri, u, v)
        cx, cy = DelaunayTriangulation.triangle_circumcenter(tri, (u, v, w))
        push!(poly_points, (cx, cy))
        v = w
    end
    vx, vy = get_point(tri, v)
    mx, my = (ux + vx) / 2, (uy + vy) / 2
    push!(poly_points, (mx, my), (mx1, my1))
    return polygon_area(poly_points)
end

function post_insertion_area(tri::Triangulation, i, envelope_vertices, point)
    u = envelope_vertices[i]
    prev_index = i == 1 ? length(envelope_vertices) - 1 : i - 1
    next_index = i == length(envelope_vertices) ? 1 : i + 1
    first_neighbour = envelope_vertices[prev_index]
    last_neighbour = envelope_vertices[next_index]
    p, q, r = get_point(tri, u, first_neighbour, last_neighbour)
    px, py = getxy(p)
    qx, qy = getxy(q)
    rx, ry = getxy(r)
    mpq = (px + qx) / 2, (py + qy) / 2
    mpr = (px + rx) / 2, (py + ry) / 2
    g1 = DelaunayTriangulation.triangle_circumcenter(p, r, point)
    !all(isfinite, g1) && return NaN # point is one of p and r, i.e. we are interpolating at a data site
    g2 = DelaunayTriangulation.triangle_circumcenter(q, p, point)
    !all(isfinite, g2) && return NaN
    points = (mpq, mpr, g1, g2, mpq)
    return polygon_area(points)
end

function compute_sibson_coordinates(tri::Triangulation, envelope_vertices, point)
    coordinates = zeros(length(envelope_vertices) - 1)
    w = 0.0
    for i in firstindex(envelope_vertices):(lastindex(envelope_vertices) - 1)
        pre = max(0.0, pre_insertion_area(tri, i, envelope_vertices))
        post = max(0.0, post_insertion_area(tri, i, envelope_vertices, point))
        if isnan(post) # need to return the the vector λ = [1] since we are exactly at a data site
            return [1.0]
        end
        coordinates[i] = max(pre - post, 0.0) # take care of any precision issues
        w += coordinates[i]
    end
    coordinates ./= w
    return coordinates
end

function evaluate_sibson_interpolant(tri::Triangulation, z, point)
    envelope_vertices, _ = compute_envelope(tri, point)
    λ = compute_sibson_coordinates(tri, envelope_vertices, point)
    if length(λ) == 1
        for i in each_solid_vertex(tri)
            get_point(tri, i) == point && return z[i]
        end
    else
        itp = 0.0
        for (λ, k) in zip(λ, envelope_vertices)
            itp += λ * z[k]
        end
        return itp
    end
end

f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)
trit = triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 30, 30)
zz = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(trit)]
xx = LinRange(0.001, 0.999, 20) # handling points on the boundary requires more care than we have discussed here
yy = LinRange(0.001, 0.999, 20)
fig = Figure(fontsize = 24)
ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"y", title = "True function", titlealign = :left, width = 400, height = 400)
contourf!(ax, xx, yy, f.(xx, yy'))
ax2 = Axis(fig[1, 2], xlabel = L"x", ylabel = L"y", title = "Interpolant", titlealign = :left, width = 400, height = 400)
zi = [evaluate_sibson_interpolant(trit, zz, (xᵢ, yᵢ)) for xᵢ in xx, yᵢ in yy]
contourf!(ax2, xx, yy, zi)
resize_to_layout!(fig)
fig

This page was generated using Literate.jl.

  • 1This is more expensive than we need. In NaturalNeighbours.jl, we use the peek keyword in triangulate to avoid making any changes to the triangulation itself, and use the InsertionEventHistory to track all changes made.