diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bade8e8..a3a050b 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/Project.toml b/Project.toml index ea3903b..50eb52a 100644 --- a/Project.toml +++ b/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" diff --git a/src/ScatteredInterpolation.jl b/src/ScatteredInterpolation.jl index 86c27fa..d0ae2c0 100755 --- a/src/ScatteredInterpolation.jl +++ b/src/ScatteredInterpolation.jl @@ -1,9 +1,15 @@ 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 @@ -11,6 +17,7 @@ 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 diff --git a/src/linear.jl b/src/linear.jl new file mode 100644 index 0000000..f2c6fde --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/linear.jl b/test/linear.jl new file mode 100644 index 0000000..1990caf --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 76b70c7..f10f819 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,3 +11,4 @@ end include("rbf.jl") include("idw.jl") include("nearestNeighbor.jl") +include("linear.jl")