2D irregular interpolation

[tomasaschan]

I have worked out a scheme for linear interpolation of functions f(x,y) that I think might be fast enough to be useful (and might be possible to optimize really well with the same metaprogramming tricks we use on B-splines). It's probably not novel, but I don't know what it's called, so I'm just going to describe it.

Given lists of points (x,y,z) (maybe on the form of equal-length lists xs, ys and zs) representing a discretization of a function z = f(x,y), we can find the value of z at an arbitrary point inside the domain using the following steps:

• Initialization:
  • Create a triangular mesh on the xy-plane
• Evaluation:

  1. Find the three corners of the triangle that contains the sought point (x,y). Let's call them P, Q and R.

  2. Evaluate z as follows:

     A = cross(Q-P, R-P)
     d = dot(A, P)
     z = (d - (x * A[1] + y * A[2])) / A[3]
     

     (These relations follow from finding the equation of the plane through the three points, and then solving for z.)

I would love to have this functionality in the package, and I could take a stab at an implementation (at least as a POC) but I have no idea how to build the triangular mesh. Maybe a Delaunay triangulation is what we want? Is there a package that does that?

[vchuravy]

There is https://github.com/JuliaGeometry/VoronoiDelaunay.jl

[tomasaschan]

@vchuravy: Nice, thanks for locating that so quickly for me! Seems like it should be able to do what we need.

[timholy]

This doesn't have to be limited to 2d: the same strategy works for polyhedra in N dimensions.

[tomasaschan]

@timholy I was hoping you'd jump onto this thread with info like that :) Do you know if this method has a name, so I can read up on it? If not, could you outline for me how to generalize the calculation given a set of n corners?

Such a generalization would also require a general implementation for finding a hyper-tetrahedron mesh, which is definitiely out-of-scope for this package (and, from what I can find with a few quick searches, doesn't seem to exist in Julia yet)...

[timholy]

I don't know if it has a name. I would search for terms like "piecewise linear polyhedra voronoi" or something. I'm attaching an old (>10 years) document I once wrote up before I realized this stuff must be well known. It was really targeted at defining PDEs on irregular grids, so it may not be entirely relevant, but for what it's worth...
delaunops.pdf

As far as finding the mesh goes, perhaps https://github.com/davidavdav/CHull.jl? (Matlab uses the QHull library for its delaunayn function, but I haven't checked wethere that's part of CHull.) CGAL also has an implementation, but that may be harder to wrap since it's C++. Given that VoronoiDelaunay claims to be faster than CGAL, obviously the best would be for someone who needs that functionality to do a julia implementation.

[yakir12]

OOO, how's this going? Can I help? I need it.. See here.

[timholy]

I'm pretty sure no one is working on this. Would be great to have someone contribute it!

[aterenin]

+1 to please add this - highly useful!

[karajan9]

I don't really have any knowledge of this (except what you can read on Wikipedia) or the internals of Interpolations, but I ported this script (which implements polyharmonic splines and works in N dims) to 1.1 and fixed the performance in the obvious places.

It's not particularly fast but it seems to work for my problems. If this is something that could fit into Interpolations I'm happy to help to make that happen -- usually need this more than the gridded alternatives.

Here is the code. Example usage would be:

f(x, y, z) = (x - y - 1.0) + exp(z)
x = rand(100)
y = rand(100)
z = rand(100)
a = f.(x, y, z)

p = PolyharmonicSpline(3, [x y z], a)  # 3 is the order

x2 = rand(10)
y2 = rand(10)
z2 = rand(10)
interpolate(p, [x2 y2 z2])

[zhiyuanzhai]

Is there any update on this?

[yakir12]

Thanks to the excellent explanation and breakdown by @tomasaschan at the beginning of this issue, this ~6 years-old issue seems astonishingly simple to solve for linear interpolations. I must be missing something or does this look good as a starting point for a PR?

using Statistics, LinearAlgebra

using Delaunay, GeometryBasics # we need a Delaunay library and something that can detect if a point is in a triangle
import GeometryBasics.Ngon # just for type parameters' sake

# this struct keeps both the 2D and 3D versions of the triangles
struct Ungridded2D{T<:Real}
    triangles::Vector{Pair{Ngon{2, T, 3, Point2{T}}, Ngon{3, T, 3, Point3{T}}}}
    function Ungridded2D{T}(xy, z) where T <: Real
        mat = hcat(first.(xy), last.(xy))
        m = delaunay(mat)
        triangles = [Triangle(xy[i]...) => Triangle((Point3{T}(xy[j]..., z[j]) for j in i)...) for i in eachrow(m.simplices)]
        return new(triangles)
    end
end
Ungridded2D(xy::Vector{Point2{T}}, z) where {T<:Real} = Ungridded2D{T}(xy, z)

function (ug::Ungridded2D)(xy)
    for (tri2, tri3) in ug.triangles
        if xy in tri2 # check if the point is inside the 2D triangle
            P, Q, R = tri3 # use the 3D triangle in @tomasaschan's algorithm
            A = (Q-P) × (R-P) 
            d = A ⋅ P
            x, y = xy
            z = (d - (x * A[1] + y * A[2])) / A[3]
            return z
        end
    end
    return NaN # defaulting to NaNs for extrapolation, for now
end

fun(xy) = 3*xy[1] - 2*xy[2] + 10 # a plane
n = 10
xy = rand(Point{2, Float64}, n)
append!(xy, [Point2(0,0), Point2(0,1), Point2(1,0), Point2(1,1)]) # required to avoid the extrapolated NaNs, just for this MWE
z = fun.(xy)
itp = Ungridded2D(xy, z)

all(1:1000) do _
    xy = rand(Point2{Float64})
    z = fun(xy)
    ẑ = itp(xy)
    isapprox(z, ẑ)
end # is true

using BenchmarkTools

@benchmark Ungridded2D($xy, $z)
# BenchmarkTools.Trial: 4895 samples with 1 evaluation.
#  Range (min … max):  974.792 μs …  16.887 ms  ┊ GC (min … max): 0.00% … 40.20%
#  Time  (median):     995.583 μs               ┊ GC (median):    0.00%
#  Time  (mean ± σ):     1.020 ms ± 498.802 μs  ┊ GC (mean ± σ):  0.66% ±  1.28%
#
#      ▂▅▆▇███▇▆▅▅▄▃▃▂▁▁▁▁  ▂▂▃▃▃▂▁▁                              ▂
#   ▃▄▆█████████████████████████████▇████▇▇▇▆▆▇▇▆▇▆▇▇▇▅▆▅▅▅▅▁▅▅▄▅ █
#   975 μs        Histogram: log(frequency) by time       1.11 ms <
#
#  Memory estimate: 47.78 KiB, allocs estimate: 1064.

@benchmark itp($(rand(Point2{Float64})))
# BenchmarkTools.Trial: 10000 samples with 976 evaluations.
#  Range (min … max):  69.117 ns …  1.763 μs  ┊ GC (min … max): 0.00% … 94.65%
#  Time  (median):     69.758 ns              ┊ GC (median):    0.00%
#  Time  (mean ± σ):   72.308 ns ± 52.368 ns  ┊ GC (mean ± σ):  2.26% ±  2.99%
#
#   ▆▆█▇▄▄▂▁ ▃▆▆▃▂ ▃▃▁    ▁▂▁     ▁▁▁ ▄▅▁                       ▂
#   █████████████████████▇███████████████▇▆▆▆▅▄▆▅▄▄▃▄▅▅▅▄▅▅▆▄▄▂ █
#   69.1 ns      Histogram: log(frequency) by time      77.3 ns <
#
#  Memory estimate: 48 bytes, allocs estimate: 2.

[mkitti]

The origin of this precedes my time as maintainer.

What I am leaning towards here is creating subpackages to avoid the dependencies of Interpolations.jl incrementally increasing over time.

Another thought is whether this better fits ScatteredInterpolations.jl.

Overall the idea for expanding Interpolations.jl would be to maintain a unified interface for multiple interpolation schemes.

In particular, we would the new interpolation type to derive from AbstractInterpolation and be able to use the other BSpline schemes in the future. If that sounds fine to you, open up a pull request and then we can discuss the details.

[mkitti]

For the subpackage, the idea would be to create a package in a lib folder. My suggestion would be to be as specific as possible when naming the package and types. In this, I would suggest including Delaunay in the name. There may be other schemes to interpolate irregular data.

[yakir12]

Hi @mkitti!
I'm not sure what the point of hosting something like this in a submodule will be (e.g. no benefit to precompilation times). I agree that it might be easier to add to ScatteredInterpolations.jl (thanks for mentioning that, I didn't know about it), if at all.

As to the original issue here, it seems like the following packages:

• https://github.com/eljungsk/ScatteredInterpolation.jl
• https://github.com/gher-ulg/DIVAnd.jl
• https://github.com/JuliaEarth/GeoStats.jl
• https://github.com/GenericMappingTools/GMT.jl
• https://github.com/kbarbary/Dierckx.jl
• (python) https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CloughTocher2DInterpolator.html

implement some kind of interpolation for un-gridded data. If Interpolations.jl was to offer 2D irregular interpolation (i.e. this issue), then Interpolations.jl would need to either depend on them or reimplement the needed functionalities in an idiomatic way to Interpolations.jl. I don't think there's a way to escape the cost of adding this feature (i.e. by keeping it in a submodule).

[mkitti]

I explained the lib folder and subpackages in this post:

https://discourse.julialang.org/t/can-someone-please-explain-the-lib-directory/84238?u=mkitti

Essentially this would be another package hosted in this repository to aid co-development.

Because it would be a separate package, it would have fully separate dependencies. It would also have independent precompilation times. However, we could fully integrate CI and the documentation.

https://discourse.julialang.org/t/usage-of-subdirectories-to-store-multiple-packages-in-a-single-repo/55534

[mkitti]

I explained the lib folder and subpackages in this post:

https://discourse.julialang.org/t/can-someone-please-explain-the-lib-directory/84238?u=mkitti

Essentially this would be another package hosted in this repository to aid co-development.

Because it would be a separate package, it would have fully separate dependencies. However, we could fully integrate CI and the documentation.

https://discourse.julialang.org/t/usage-of-subdirectories-to-store-multiple-packages-in-a-single-repo/55534

[mkitti]

A future direction of Interpolations.jl might be to break it up into smaller packages, and provide a meta package that depends on all of them. Perhaps Interpolations.jl becomes the meta package in the future.

[fatteneder]

(python) https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CloughTocher2DInterpolator.html

We recently ported this one to Julia https://github.com/fatteneder/CloughTocher2DInterpolation.jl and will soon be using it in https://github.com/MakieOrg/TopoPlots.jl.