Laplace

src/NaturalNeighbourInterp.jl

@@ -65,6 +65,7 @@ include("coordinates/main.jl")
 include("coordinates/sibson.jl")
 include("coordinates/triangle.jl")
 include("coordinates/nearest.jl")
+include("coordinates/laplace.jl")
 include("coordinates/extrapolation.jl")
 include("interpolate.jl")
 include("utils.jl")

src/coordinates/laplace.jl

@@ -0,0 +1,43 @@
+function _compute_laplace_coordinates(
+    tri::Triangulation{P,Ts,I,E,Es,BN,BNM,B,BIR,BPL},
+    interpolation_point,
+    cache::InterpolantCache{F}=InterpolantCache(tri);
+    kwargs...
+) where {P,Ts,I,E,Es,BN,BNM,B,BIR,BPL,F}
+    coordinates = get_coordinates(cache)
+    envelope = get_envelope(cache)
+    insertion_event_history = get_insertion_event_history(cache)
+    temp_adjacent = get_temp_adjacent(cache)
+    last_triangle = get_last_triangle(cache)
+    envelope, temp_adjacent, insertion_event_history, V = compute_bowyer_envelope!(envelope, tri, insertion_event_history, temp_adjacent, interpolation_point; try_points=last_triangle[], kwargs...) #kwargs are add_point! kwargs
+    i, j, return_flag = check_for_extrapolation(tri, V, interpolation_point, last_triangle)
+    return_flag && return two_point_interpolate!(coordinates, envelope, tri, i, j, interpolation_point)
+    resize!(coordinates, length(envelope) - 1)
+    w = zero(number_type(tri))
+    for i in firstindex(envelope):(lastindex(envelope)-1)
+        ratio = laplace_ratio(tri, envelope, i, interpolation_point) # could reuse a circumcenter here, but it's not the dominating part of the computation anyway.
+        isnan(ratio) && return handle_duplicate_points!(tri, interpolation_point, coordinates, envelope)
+        coordinates[i] = ratio
+        w += coordinates[i]
+    end
+    pop!(envelope)
+    coordinates ./= w
+    return NaturalCoordinates(coordinates, envelope, interpolation_point, tri)
+end
+
+function laplace_ratio(tri, envelope, i, interpolation_point)
+    u = envelope[i]
+    prev_u = envelope[previndex_circular(envelope, i)]
+    next_u = envelope[nextindex_circular(envelope, i)]
+    p, q, r = get_point(tri, u, prev_u, next_u)
+    g1 = triangle_circumcenter(q, p, interpolation_point)
+    g2 = triangle_circumcenter(p, r, interpolation_point)
+    g1x, g1y = getxy(g1)
+    g2x, g2y = getxy(g2)
+    ℓ² = (g1x - g2x)^2 + (g1y - g2y)^2
+    px, py = getxy(p)
+    x, y = getxy(interpolation_point)
+    d² = (px - x)^2 + (py - y)^2
+    w = sqrt(ℓ² / d²)
+    return w
+end

src/coordinates/main.jl

@@ -11,7 +11,9 @@ function compute_natural_coordinates(
         return _compute_triangle_coordinates(tri, interpolation_point, cache; kwargs...)
     elseif method == :nearest # not local coordinates, but still a nice method to have...
         return _compute_nearest_coordinates(tri, interpolation_point, cache; kwargs...)
+    elseif method == :laplace
+        return _compute_laplace_coordinates(tri, interpolation_point, cache; kwargs...)
     else
-        throw(ArgumentError("method must be one of :sibson, :triangle, or :nearest."))
+        throw(ArgumentError("method must be one of :sibson, :triangle, :laplace, or :nearest."))
     end
 end

src/coordinates/sibson.jl

@@ -15,25 +15,15 @@ function _compute_sibson_coordinates(
     return_flag && return two_point_interpolate!(coordinates, envelope, tri, i, j, interpolation_point)
     resize!(coordinates, length(envelope) - 1)
     w = zero(number_type(tri))
-    duplicate_point_flag = false
     for i in firstindex(envelope):(lastindex(envelope)-1)
         pre = pre_insertion_area!(poly_points, envelope, i, tri)
         post = post_insertion_area(envelope, i, tri, interpolation_point)
-        if isnan(post)
-            duplicate_point_flag = true
-            break
-        end
+        isnan(post) && return handle_duplicate_points!(tri, interpolation_point, coordinates, envelope)
         coordinates[i] = pre - post
         w += coordinates[i]
     end
     pop!(envelope)
-    if duplicate_point_flag
-        idx = findfirst(i -> get_point(tri, i) == getxy(interpolation_point), envelope)
-        fill!(coordinates, zero(F))
-        coordinates[idx] = one(F)
-    else
-        coordinates ./= w
-    end
+    coordinates ./= w
     return NaturalCoordinates(coordinates, envelope, interpolation_point, tri)
 end
 

src/coordinates/triangle.jl

@@ -7,7 +7,7 @@ function _compute_triangle_coordinates(
     coordinates = get_coordinates(cache)
     envelope = get_envelope(cache)
     last_triangle = get_last_triangle(cache)
-    V = jump_and_march(tri, interpolation_point; try_points=last_triangle[])
+    V = jump_and_march(tri, interpolation_point; try_points=last_triangle[], kwargs...)
     i, j, return_flag = check_for_extrapolation(tri, V, interpolation_point, last_triangle)
     return_flag && return two_point_interpolate!(coordinates, envelope, tri, i, j, interpolation_point)
     i, j, k = indices(V)

src/utils.jl

@@ -66,4 +66,14 @@ function get_barycentric_deviation(natural_coordinates::NaturalCoordinates{F}) w
     x, y = getxy(interpolation_point)
     δ² = (x - x̂)^2 + (y - ŷ)^2
     return sqrt(δ²)
+end
+
+function handle_duplicate_points!(tri, interpolation_point, coordinates::AbstractVector{F}, envelope) where {F}
+    idx = findfirst(i -> get_point(tri, i) == getxy(interpolation_point), envelope)
+    envelope_idx = envelope[idx]
+    resize!(coordinates, 1)
+    resize!(envelope, 1)
+    envelope[begin] = envelope_idx
+    coordinates[begin] = one(F)
+    return NaturalCoordinates(coordinates, envelope, interpolation_point, tri)
 end

test/runtests.jl

@@ -283,7 +283,7 @@ end
 end
 
 @testset "Natural coordinates" begin
-    for method in (:sibson, :triangle, :nearest)
+    for method in (:sibson, :triangle, :nearest, :laplace)
         pts = [(0.0, 8.0), (0.0, 0.0), (14.0, 0.0), (14.0, 8.0), (4.0, 4.0), (10.0, 6.0), (6.0, 2.0), (12.0, 4.0), (0.0, 4.0)]
         tri = triangulate(pts, randomise=false, delete_ghosts=false)
         n = 2500
@@ -344,7 +344,7 @@ end
 end
 
 function test_interpolant(itp, x, y, f)
-    for method in (:sibson, :triangle, :nearest)
+    for method in (:sibson, :triangle, :nearest, :laplace)
         for _ in 1:500
             vals = itp(x, y; parallel=false, method)
             vals2 = similar(vals)
@@ -421,9 +421,12 @@ end
     t = dbp / dab
     _z = t * z[i] + (1 - t) * z[j]
     __z = itp(getx(p), gety(p); method=:triangle)
-    @test _z ≈ __z
+    @test _z ≈ itp(getx(p), gety(p); method=:triangle)
+    @test _z ≈ itp(getx(p), gety(p); method=:sibson)
+    @test _z ≈ itp(getx(p), gety(p); method=:laplace)
     @test __z ≈ itp(1.8, 0.7; method=:triangle)
     @test __z ≈ itp(1.8, 0.7; method=:sibson)
+    @test __z ≈ itp(1.8, 0.7; method=:laplace)
 end
 
 @testset "Example I: No extrapolation" begin
@@ -536,4 +539,76 @@ end
     end
     resize_to_layout!(fig)
     @test_reference "figures/example_2.png" fig
+end
+
+@testset "Test coefficient values for each method" begin # used GeoGebra
+    for _ in 1:100 # make sure rng is getting passed consistently 
+        # Setup
+        rng = StableRNG(872973)
+        pts = [(0.0, 8.0), (0.0, 0.0), (14.0, 0.0),
+            (14.0, 8.0), (4.0, 4.0), (10.0, 6.0),
+            (6.0, 2.0), (12.0, 4.0), (0.0, 4.0),
+            (2.5, 5.0), (7.0, 3.3), (4.5, 5.2),
+            (13.0, 0.5), (12.0, 6.0), (8.5, 3.5),
+            (0.5, 6.0), (1.5, 6.0), (3.5, 6.0),
+            (0.5, 2.0), (2.5, 2.0), (2.5, 2.5),
+            (9.0, 2.0), (8.5, 6.0), (4.0, 2.0)]
+        tri = triangulate(pts, randomise=false, delete_ghosts=false, rng=rng)
+        vorn = voronoi(tri, false)
+        q = (5.0, 4.0)
+        tri2 = deepcopy(tri)
+        add_point!(tri2, q, rng=rng)
+        vorn2 = voronoi(tri2, false)
+        V = get_polygon(vorn2, num_points(tri2))
+        AX2 = get_area(vorn2, num_points(tri2))
+
+        # Sibson
+        nc = NNI.compute_natural_coordinates(tri, q; method=:sibson, rng=rng)
+        AF1G1D1B1A1 = 1.1739978952813 # E = 5
+        A1B1C1 = 0.062500000375 # Z = 24
+        B1D1E1C1 = 0.2540749084958 # G = 7
+        H1I1E1D1G1 = 0.7376579777378 # K = 11
+        F1G1J1K1 = 0.9489911839831 # L = 12
+        K1I1J1 = 0.003313286868 # W = 23
+        AX = AF1G1D1B1A1 + A1B1C1 + B1D1E1C1 + H1I1E1D1G1 + F1G1J1K1 + K1I1J1
+        @test AX ≈ AX2 rtol = 1e-3
+        @test nc.indices == [23, 12, 5, 24, 7, 11]
+        @test nc.coordinates ≈ [K1I1J1, F1G1J1K1, AF1G1D1B1A1, A1B1C1, B1D1E1C1, H1I1E1D1G1] ./ AX rtol = 1e-2
+
+        # Triangle 
+        nc = NNI.compute_natural_coordinates(tri, q; method=:triangle, rng=rng)
+        V = jump_and_march(tri, q; rng)
+        @test nc.indices == [5, 11, 12]
+        @test nc.coordinates ≈ [0.52, 0.3, 0.18] rtol = 1e-2
+
+        # Laplace 
+        nc = NNI.compute_natural_coordinates(tri, q; method=:laplace, rng=rng)
+        dqw = 4.0311288741493
+        dqk = 2.1189620100417
+        dqg = 2.2360679774998
+        dqz = 2.2360679774998
+        dqe = 1.0
+        dqℓ = 1.3
+        dc1b1 = 2.2893491697301
+        da1b1 = 0.0572608008105
+        df1b1 = 2.2260773066834
+        df1e1 = 1.4888476232694
+        de1d1 = 0.5650898843856
+        dd1c1 = 0.9335156761474
+        k = dc1b1 / dqk
+        w = da1b1 / dqw
+        ℓ = df1b1 / dqℓ
+        e = df1e1 / dqe
+        z = de1d1 / dqz
+        g = dd1c1 / dqg
+        tot = k + w + ℓ + e + z + g
+        k /= tot
+        w /= tot
+        ℓ /= tot
+        e /= tot
+        z /= tot
+        g /= tot
+        @test nc.indices == [23, 12, 5, 24, 7, 11]
+        @test nc.coordinates ≈ [w, ℓ, e, z, g, k] rtol = 1e-2
+    end
 end