From 4cb4553c94294d4fdb5c2bee6a469d1506741dad Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Sun, 8 Aug 2021 08:40:13 +0200 Subject: [PATCH 1/7] Initial implementation of linear interpolation --- Project.toml | 6 +++ src/ScatteredInterpolation.jl | 13 ++++-- src/linear.jl | 81 +++++++++++++++++++++++++++++++++++ 3 files changed, 97 insertions(+), 3 deletions(-) create mode 100644 src/linear.jl diff --git a/Project.toml b/Project.toml index ea3903b..09159c3 100644 --- a/Project.toml +++ b/Project.toml @@ -3,15 +3,21 @@ 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" +StaticArrays = "1" julia = "1" [extras] 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..e301fc4 --- /dev/null +++ b/src/linear.jl @@ -0,0 +1,81 @@ + +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 + @assert eltype(points) == Float64 "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 + 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) + + for (i, bs) in enumerate(itp.barycentric_setup) + coord = cartesian2barycentric(bs, point) + + if all(coord .≥ 0) && all(coord .≤ 1) + return (i, coord) + end + end + + error("Data out of range at $point. Extrapolation is not supported.") +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 From 0979f2ba9d84ae7f500e04ac7ce19ebafbd51f93 Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Sun, 8 Aug 2021 08:53:05 +0200 Subject: [PATCH 2/7] Restrict minimum Julia version to 1.4 (requirement by Delaunay.jl) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 09159c3..50eb52a 100644 --- a/Project.toml +++ b/Project.toml @@ -18,7 +18,7 @@ Delaunay = "1" Distances = "0.9,0.10" NearestNeighbors = "0.4" StaticArrays = "1" -julia = "1" +julia = "1.4" [extras] Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" From b23f0a856b5c0c5cc9e18fdb5fe5dbfddfc820bd Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Sun, 8 Aug 2021 08:54:48 +0200 Subject: [PATCH 3/7] Stop running CI on 1.0 --- .github/workflows/ci.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bade8e8..07fcb2c 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: From dd84d99525122b9adeed70544a01b2ae07cb7437 Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Sun, 8 Aug 2021 09:00:55 +0200 Subject: [PATCH 4/7] Update ci.yml --- .github/workflows/ci.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 07fcb2c..a3a050b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -38,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 @@ -51,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 From 75b6e07b3d0265ad0d7c54429621dcffb31f2434 Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Sun, 10 Oct 2021 07:33:39 +0200 Subject: [PATCH 5/7] Update error handling --- src/linear.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/linear.jl b/src/linear.jl index e301fc4..379c55c 100644 --- a/src/linear.jl +++ b/src/linear.jl @@ -19,7 +19,7 @@ function interpolate(m::Linear, samples::AbstractArray{<:Number, N}) where {N} # Delaunay.jl requires Float64 for the points at the moment - @assert eltype(points) == Float64 "Linear interpolation is only supported for Float64 point coordinates" + 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)) @@ -46,6 +46,9 @@ function evaluate(itp::LinearInterpolant, points::AbstractArray{<:Real, 2}) 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) @@ -71,7 +74,7 @@ function find_simplex(itp, point) end end - error("Data out of range at $point. Extrapolation is not supported.") + throw(DomainError(point, "Extrapolation is not supported.")) end # Helper function to copy the points array into a vector of SVector. From 4360df00d893cc79a67807dfb042ec9db4b18b20 Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Sun, 10 Oct 2021 07:53:51 +0200 Subject: [PATCH 6/7] Add tests --- test/linear.jl | 25 +++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 26 insertions(+) create mode 100644 test/linear.jl 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") From d71dfb6b6276bbaf252c9f9aefed780c71267abc Mon Sep 17 00:00:00 2001 From: Emil Ljungskog Date: Mon, 18 Jul 2022 17:29:18 +0200 Subject: [PATCH 7/7] Allow for rounding error on domain edges. --- src/linear.jl | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/src/linear.jl b/src/linear.jl index 379c55c..f2c6fde 100644 --- a/src/linear.jl +++ b/src/linear.jl @@ -66,15 +66,33 @@ 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 .≥ 0) && all(coord .≤ 1) + if all(coord .≥ -tol) && all(coord .≤ 1 + 2*tol) return (i, coord) end end - - throw(DomainError(point, "Extrapolation is not supported.")) + return nothing end # Helper function to copy the points array into a vector of SVector.