Linear interpolation

.github/workflows/ci.yml

@@ -15,7 +15,6 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
           - '1'
           - 'nightly'
         os:
@@ -39,6 +38,8 @@ jobs:
             ${{ runner.os }}-test-
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
+        env:
+          PYTHON:
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v1
@@ -52,6 +53,9 @@ jobs:
       - uses: julia-actions/setup-julia@v1
         with:
           version: '1'
+      - uses: julia-actions/julia-buildpkg@v1
+        env:
+          PYTHON:
       - run: |
           julia --project=docs -e '
             using Pkg

Project.toml

@@ -3,16 +3,22 @@ uuid = "3f865c0f-6dca-5f4d-999b-29fe1e7e3c92"
 version = "0.3.6"
 
 [deps]
+Bernstein = "e9b0fb4c-9cb7-4f61-9c14-701a41c684d7"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+Delaunay = "07eb4e4e-0c6d-46ef-bc4e-83d5e5d860a9"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+Bernstein = "1"
 Combinatorics = "1"
+Delaunay = "1"
 Distances = "0.9,0.10"
 NearestNeighbors = "0.4"
-julia = "1"
+StaticArrays = "1"
+julia = "1.4"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

src/ScatteredInterpolation.jl

@@ -1,16 +1,23 @@
 module ScatteredInterpolation
 
-using Distances, NearestNeighbors, Combinatorics, LinearAlgebra
+using   Distances,
+        NearestNeighbors,
+        Combinatorics,
+        LinearAlgebra,
+        Delaunay,
+        StaticArrays,
+        Bernstein
 
-export interpolate,
-    evaluate
+export  interpolate,
+        evaluate
 
 abstract type ScatteredInterpolant end
 abstract type InterpolationMethod end
 
 include("./rbf.jl")
 include("./idw.jl")
 include("./nearestNeighbor.jl")
+include("./linear.jl")
 
 
 # Fallback method for the case of just one point

src/linear.jl

@@ -0,0 +1,102 @@
+
+export Linear
+
+"""
+    Linear()
+
+Linear interpolation based on the Delaunay triangulation of the sample points.
+
+"""
+struct Linear <: InterpolationMethod end
+
+struct LinearInterpolant{S, BT} <: ScatteredInterpolant where {S <: AbstractArray{<:Number, N}, BT} where {N}
+    samples::S
+    barycentric_setup::BT
+end
+
+function interpolate(m::Linear,
+            points::AbstractArray{<:Real, 2},
+            samples::AbstractArray{<:Number, N}) where {N}
+
+    # Delaunay.jl requires Float64 for the points at the moment
+    eltype(points) == Float64 || throw(ArgumentError("Linear interpolation is only supported for Float64 point coordinates"))
+
+    dim = size(points, 1)
+    spoints = copy_svec(Float64, points, Val(dim))
+
+    # Build the Delaunay triangulation
+    # Delaunay.jl expects the points array to be nPoints x nDims, we take the opposite as input
+    points = permutedims(points)
+    mesh = delaunay(points)
+
+    # Extract each simplex as a coordinate matrix with one point per row, and prepare for computing the Barycentric coordinates
+    nSimplices = size(mesh.simplices, 1)
+    simplices = [SMatrix{dim+1, dim, eltype(points)}(points[mesh.simplices[i, :], :]) for i in 1:nSimplices]
+    barycentric_setup = cartesian2barycentric_setup.(simplices)
+
+    # Group the samples for each simplex
+    grouped_samples = [SMatrix{dim+1, size(samples, 2), eltype(samples)}(samples[mesh.simplices[i, :], :]) for i in 1:nSimplices]
+
+    return LinearInterpolant(grouped_samples, barycentric_setup)
+end
+
+function evaluate(itp::LinearInterpolant, points::AbstractArray{<:Real, 2})
+
+    # Compute a reasonable type for the output data
+    T = StaticArrays.arithmetic_closure(eltype(itp.samples[1]))
+
+    dim = size(itp.samples[1], 1) - 1
+    pointsDim = size(points, 1)
+    dim == pointsDim || throw(DimensionMismatch("got points in dimension $pointsDim but expected $dim."))
+
+    spoints = copy_svec(eltype(points), points, Val(dim))
+
+    m = length(spoints)
+    n = size(itp.samples[1], 2)
+    values = zeros(T, m, n)
+
+    for (i, point) in enumerate(spoints)
+        (ind, bCoord) = find_simplex(itp, point)
+        values[i,:] = bCoord'*itp.samples[ind]
+    end
+
+    return values
+end
+
+# Find the simplex containing the point to interpolate, and return the barycentric coordinates
+function find_simplex(itp, point)
+
+    coord = simplex_search(itp, point, 0.0)
+    if coord !== nothing
+        return coord
+    end
+
+    # Look again, but this time with a small tolerance to allow for rounding errors at the
+    # edge of the domain.
+    coord = simplex_search(itp, point, 5*eps())
+    if coord !== nothing
+        return coord
+    end
+
+    throw(DomainError(point, "Extrapolation is not supported."))
+end
+
+# Find the simplex containing the point to interpolate, and return the barycentric coordinates
+# Use a tolerance to allow for some rounding error at the edge of the domain.
+@inline function simplex_search(itp, point, tol)
+
+    for (i, bs) in enumerate(itp.barycentric_setup)
+        coord = cartesian2barycentric(bs, point)
+
+        if all(coord .≥ -tol) && all(coord .≤ 1 + 2*tol)
+            return (i, coord)
+        end
+    end
+    return nothing
+end
+
+# Helper function to copy the points array into a vector of SVector.
+# Borrowed from NearestNeighbors.jl
+@inline function copy_svec(::Type{T}, data, ::Val{dim}) where {T, dim} 
+    [SVector{dim, T}(ntuple(i -> data[n+i], Val(dim))) for n in 0:dim:(length(data)-1)]
+end

test/linear.jl

@@ -0,0 +1,25 @@
+# Define some points and data in 2D
+points = permutedims([0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0], (2,1))
+points_adjoint = [0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0]'
+data = [0.0; 0.5; 0.5; 1.0]
+
+@testset "Linear" begin
+
+    @testset "Evaluation" for p in (points, points_adjoint)
+
+        itp = interpolate(Linear(), p, data)
+
+        # Check that we get back the original data at the sample points
+        ev = evaluate(itp, points)
+        @test ev ≈ data
+
+        # Make sure that we cannot do extrapolation
+        @test_throws DomainError evaluate(itp, [2.0; 2.0]) 
+        @test_throws DomainError evaluate(itp, [-1.0; -1.0]) 
+
+        # Test for linearity
+        tp = [0.2 0.2; 0.4 0.4; 0.6 0.6; 0.8 0.8]'
+        ev = evaluate(itp, tp)
+        @test ev ≈ tp[1,:]
+    end
+end

test/runtests.jl

@@ -11,3 +11,4 @@ end
 include("rbf.jl")
 include("idw.jl")
 include("nearestNeighbor.jl")
+include("linear.jl")