Merge pull request #19 from DanielVandH/allow-cons

Allow for constrained triangulations

Project.toml

@@ -1,7 +1,7 @@
 name = "NaturalNeighbours"
 uuid = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

src/NaturalNeighbours.jl

@@ -23,6 +23,7 @@ import DelaunayTriangulation: DelaunayTriangulation,
     previndex_circular,
     nextindex_circular,
     get_point,
+    get_points,
     triangle_circumcenter,
     num_points,
     number_type,
@@ -54,7 +55,8 @@ import DelaunayTriangulation: DelaunayTriangulation,
     jump_to_voronoi_polygon,
     iterated_neighbourhood,
     iterated_neighbourhood!,
-    triangle_area
+    triangle_area,
+    get_boundary_nodes
 import ChunkSplitters: chunks
 using ElasticArrays
 using LinearAlgebra

src/data_structures/interpolant.jl

@@ -16,7 +16,7 @@ struct NaturalNeighboursInterpolant{T<:Triangulation,F,G,H,N,D}
         @assert num_solid_vertices(tri) == length(z) "The number of points in the triangulation must equal the length of the data vector."
         !has_ghost_triangles(tri) && add_ghost_triangles!(tri)
         if has_boundary_nodes(tri)
-            throw(ArgumentError("Natural neighbour interpolation is only defined over unconstrained triangulations."))
+            @warn "Natural neighbour interpolation is only defined over unconstrained triangulations. You may find unexpected results when interpolating near the boundaries or constrained edges, and especially near non-convex boundaries or outside of the triangulation. In your later evaluations of this interpolant, consider using project=false."
         end
         nt = Base.Threads.nthreads()
         derivative_caches = [DerivativeCache(tri) for _ in 1:nt]

src/interpolation/interpolate.jl

@@ -27,7 +27,7 @@ In each method, `method` defines the method used for evaluating the interpolant,
 method, `parallel` is ignored, but for the latter two methods it defines whether to use multithreading or not for evaluating the interpolant at 
 all the points. The `kwargs...` argument is passed into `add_point!` from DelaunayTriangulation.jl, e.g. you could pass some `rng`. Lastly, 
 the `project` argument determines whether extrapolation is performed by projecting any exterior points onto the boundary of the convex hull 
-of the data sites and performing two-point interpolation, or to simply replaced any extrapolated values with `Inf`.
+of the data sites and performing two-point interpolation, or to simply replace any extrapolated values with `Inf`.
 
 !!! warning
 

src/interpolation/utils.jl

@@ -6,6 +6,9 @@ function compute_bowyer_envelope!(envelope, tri::Triangulation, history::Inserti
     I = integer_type(tri)
     V = add_point!(tri, point; store_event_history=Val(true), event_history=history, peek=Val(true), kwargs...)
     all_triangles = each_added_triangle(history)
+    if isempty(all_triangles) # This is possible for constrained triangulations
+        return envelope, temp_adjacent, history, V
+    end
     for T in all_triangles
         add_triangle!(temp_adjacent, T)
     end

src/utils.jl

@@ -4,6 +4,8 @@
 Given a polygon described by `(points, boundary_nodes)`, matching the 
 specification of polygons in DelaunayTriangulation.jl (see [here](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/boundary_handling/)),
 returns a vector of indices of the points defined by `(x, y)` that are outside of the polygon.
+
+Use `tol` to specify a tolerance for the distance to the polygon.
 """
 function identify_exterior_points(x, y, points, boundary_nodes; tol = 0.0)
     @assert length(x) == length(y) "x and y must have the same length."
@@ -19,4 +21,23 @@ function identify_exterior_points(x, y, points, boundary_nodes; tol = 0.0)
         end
     end
     return exterior_points
+end
+
+"""
+    identify_exterior_points(x, y, itp::NaturalNeighboursInterpolant; tol = 0.0)
+
+Returns the indices of the points defined by the vectors `(x, y)` that are 
+outside of the underlying triangulation to the interpolant `itp`.
+
+Use `tol` to specify a tolerance for the distance to the triangulation.
+"""
+function identify_exterior_points(x, y, itp::NaturalNeighboursInterpolant; tol=0.0)
+    tri = get_triangulation(itp)
+    points = get_points(tri)
+    if !has_boundary_nodes(tri)
+        bn = get_convex_hull_indices(tri)
+    else 
+        bn = get_boundary_nodes(tri)
+    end
+    return identify_exterior_points(x, y, points, bn; tol=tol)
 end

test/interpolation/basic_tests.jl

@@ -37,9 +37,6 @@ include(normpath(@__DIR__, "../.", "helper_functions", "test_functions.jl"))
     @test NNI.get_gradient(_itp) === nothing
     @test NNI.get_hessian(_itp) === nothing
     @test itp isa NNI.NaturalNeighboursInterpolant
-    DT.lock_convex_hull!(tri)
-    @test_throws ArgumentError interpolate(tri, z)
-    DT.unlock_convex_hull!(tri)
     @test_throws AssertionError interpolate(tri, z[1:end-1])
     w = rand(length(z))
     y = rand(length(z))

test/interpolation/constrained.jl

@@ -0,0 +1,104 @@
+using ..NaturalNeighbours
+const NNI = NaturalNeighbours
+using DelaunayTriangulation
+const DT = DelaunayTriangulation
+using LinearAlgebra
+using ReferenceTests
+using CairoMakie
+function plot_2d(fig, i, j, title, vals, xg, yg, x, y, show_scatter=true)
+    ax = Axis(fig[i, j], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
+    contourf!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
+    show_scatter && scatter!(ax, vec([(i - 1) / 9 for i in (1, 3, 4, 5, 8, 9, 10), j in (1, 2, 3, 5, 6, 7, 9, 10)]), vec([(j - 1) / 9 for i in (1, 3, 4, 5, 8, 9, 10), j in (1, 2, 3, 5, 6, 7, 9, 10)]), color=:red, markersize=14)
+end
+function plot_3d(fig, i, j, title, vals, xg, yg)
+    ax = Axis3(fig[i, j], xlabel="x", ylabel="y", width=600, height=600, title=title, titlealign=:left)
+    surface!(ax, xg, yg, reshape(vals, (length(xg), length(yg))), color=vals, colormap=:viridis, levels=-1:0.05:0, extendlow=:auto, extendhigh=:auto)
+end
+
+@testset "A domain with no holes" begin
+    a, b, c, d = 0.0, 1.0, 0.0, 1.0
+    nx, ny = 10, 10
+    tri = triangulate_rectangle(a, b, c, d, nx, ny)
+    f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)
+    z = [f(x, y) for (x, y) in each_point(tri)]
+    itp = interpolate(tri, z; derivatives=true)
+    xg = LinRange(0, 1, 100)
+    yg = LinRange(0, 1, 100)
+    _x = vec([x for x in xg, _ in yg])
+    _y = vec([y for _ in xg, y in yg])
+    exact = f.(_x, _y)
+    sibson_vals = itp(_x, _y; method=Sibson())
+    triangle_vals = itp(_x, _y; method=Triangle())
+    laplace_vals = itp(_x, _y; method=Laplace())
+    sibson_1_vals = itp(_x, _y; method=Sibson(1))
+    nearest_vals = itp(_x, _y; method=Nearest())
+    farin_vals = itp(_x, _y; method=Farin())
+    hiyoshi_vals = itp(_x, _y; method=Hiyoshi(2))
+    all_vals = (sibson_vals, triangle_vals, laplace_vals, sibson_1_vals, nearest_vals, farin_vals, hiyoshi_vals, exact)
+    titles = ("(a): Sibson", "(b): Triangle", "(c): Laplace", "(d): Sibson-1", "(e): Nearest", "(f): Farin", "(g): Hiyoshi", "(h): Exact")
+    fig = Figure(fontsize=55)
+    for (i, (vals, title)) in enumerate(zip(all_vals, titles))
+        plot_2d(fig, 1, i, title, vals, xg, yg, first.(each_point(tri)), last.(each_point(tri)), !(vals == exact))
+        plot_3d(fig, 2, i, " ", vals, xg, yg)
+    end
+    resize_to_layout!(fig)
+    @test_reference normpath(@__DIR__, "../..", "example.png") fig by = psnr_equality(20)
+end
+
+@testset "A domain with holes" begin
+    R₁ = 2.0
+    R₂ = 3.0
+    θ = (collect ∘ LinRange)(0, 2π, 250)
+    θ[end] = θ[begin]
+    x = [
+        [R₂ .* cos.(θ)],
+        [reverse(R₁ .* cos.(θ))] # inner boundaries are clockwise
+    ]
+    y = [
+        [R₂ .* sin.(θ)],
+        [reverse(R₁ .* sin.(θ))] # inner boundaries are clockwise
+    ]
+    boundary_nodes, points = convert_boundary_points_to_indices(x, y)
+    tri = triangulate(points; boundary_nodes)
+    A = get_total_area(tri)
+    refine!(tri; max_area=1e-4A)
+    D = 6.25e-4
+    Tf = (x, y) -> let r = sqrt(x^2 + y^2)
+        (R₂^2 - r^2) / (4D) + R₁^2 * log(r / R₂) / (2D)
+    end
+    _safe_Tf = (x, y) -> let r = sqrt(x^2 + y^2)
+        !(R₁ ≤ r ≤ R₂) && return Inf
+        return Tf(x, y)
+    end
+    z = [Tf(x, y) for (x, y) in each_point(tri)]
+    x = first.(each_point(tri))
+    y = last.(each_point(tri))
+    triangles = [T[j] for T in each_solid_triangle(tri), j in 1:3]
+    itp = interpolate(tri, z; derivatives=true)
+    xg = LinRange(-R₂, R₂, 250)
+    yg = LinRange(-R₂, R₂, 250)
+    _x = vec([x for x in xg, _ in yg])
+    _y = vec([y for _ in xg, y in yg])
+    exact = _safe_Tf.(_x, _y)
+    sibson_vals = itp(_x, _y; method=Sibson(), project=false)
+    triangle_vals = itp(_x, _y; method=Triangle(), project=false)
+    laplace_vals = itp(_x, _y; method=Laplace(), project=false)
+    sibson_1_vals = itp(_x, _y; method=Sibson(1), project=false)
+    nearest_vals = itp(_x, _y; method=Nearest(), project=false)
+    farin_vals = itp(_x, _y; method=Farin(), project=false)
+    hiyoshi_vals = itp(_x, _y; method=Hiyoshi(2), project=false)
+    all_vals = (sibson_vals, triangle_vals, laplace_vals, sibson_1_vals, nearest_vals, farin_vals, hiyoshi_vals, exact)
+    titles = ("(a): Sibson", "(b): Triangle", "(c): Laplace", "(d): Sibson-1", "(e): Nearest", "(f): Farin", "(g): Hiyoshi", "(h): Exact")
+    _tri = triangulate([_x'; _y'])
+    _triangles = [T[j] for T in each_solid_triangle(_tri), j in 1:3]
+    fig = Figure(fontsize=55)
+    for (i, (vals, title)) in enumerate(zip(all_vals, titles))
+        ax = Axis(fig[1, i], width=600, height=600, title=title, titlealign=:left)
+        _vals = copy(vals)
+        _vals[isinf.(vals)] .= Inf
+        contourf!(ax, _x, _y, _vals, levels=0:50:900)
+    end
+    resize_to_layout!(fig)
+    fig
+    @test_reference normpath(@__DIR__, "example_constrained.png") fig
+end

test/runtests.jl

@@ -26,6 +26,9 @@ using SafeTestsets
     #@safetestset "Influence" begin
     #    include("interpolation/influence.jl")
     #end
+    @safetestset "Constrained Triangulations" begin
+        include("interpolation/constrained.jl")
+    end
 end
 
 @testset "Differentiation" begin
@@ -56,11 +59,11 @@ end
     @safetestset "Switzerland" begin
         include("doc_examples/swiss.jl")
     end
-    if get(ENV, "CI", "false") == "false"
-        @safetestset "Comparison" begin
-            include("doc_examples/interpolant_comparisons.jl")
-        end
-    end
+    #if get(ENV, "CI", "false") == "false"
+    #    @safetestset "Comparison" begin
+    #        include("doc_examples/interpolant_comparisons.jl")
+    #    end
+    #end
 end