Unwanted behaviour boundaries

[Kevin-Mattheus-Moerman]

Hi @DanielVandH, I absolute love this package well done!. I think it is vital for geometry processing etc and for scattered data interpolation/extrapolation.

I've tested this package and it works really well when interpolating well within the convexhull. And when I say "well within" I suppose that means for regions with a lot of points surrounding them rather than at the boundary.

Here is a screenshot of the performance of the various methods where I test it very much like your demo in the README but with non-gridded (random) data point input:

It looks like at the boundary we snap to a kind of "nearest neighbour/linear extrapolation origami". There are clear creases and kinks and it seems that this doesn't just happen outside of the convex hull but also just at the boundary (see the kinks starting inside the convex hull e.g. at the top corner). I suspect there are two things going on here. 1) you probably switch method for extrapolation, and 2) perhaps there are to few points at the boundary to enforce the smoothness. Would it help to reflect points out at the boundary so the local slopes continue and additional fictitious data is provided. Or is there another way, especially for the smoother methods, to make the behaviour at the boundaries more consistent with the otherwise C1/ C2 continuous nature of them?

naturalNeighbourInterp.mp4

Thanks!

For completeness, here is the code I used:

using GLMakie
using GLMakie.GeometryBasics
using NaturalNeighbours
using StableRNGs 

f = (x, y) -> sin(x * y) - cos(x - y) * exp(-(x - y)^2)

numPointsData = 100
rng = StableRNG(123)
x = rand(rng, numPointsData)
y = rand(rng, numPointsData)
z = f.(x,y)

itp = interpolate(x,y, z; derivatives=true)

numPointsInterp = 100
x_i_r = range(0.0, 1.0, numPointsInterp)
y_i_r = range(0.0, 1.0, numPointsInterp)

function pointGrid(x_i_r, y_i_r)
    x_i = Vector{Float64}(undef,numPointsInterp^2)
    y_i = Vector{Float64}(undef,numPointsInterp^2)
    i = 1
    for y in y_i_r
        for x in x_i_r        
            x_i[i] = x
            y_i[i] = y
            i+=1
        end
    end
    return x_i, y_i
end

x_i, y_i = pointGrid(x_i_r, y_i_r)

methodSet = (NaturalNeighbours.Sibson(), 
             NaturalNeighbours.Triangle(),
             NaturalNeighbours.Laplace(),
             NaturalNeighbours.Sibson(1),
             NaturalNeighbours.Nearest(),
             NaturalNeighbours.Farin(),
             NaturalNeighbours.Hiyoshi(1),
             NaturalNeighbours.Hiyoshi(2)
             )

titleStrings = ("Sibson", "Triangle", "Laplace", "Sibson(1)", "Nearest", "Farin", "Hiyoshi(1)", "Hiyoshi(2)")             
fig = Figure(size = (800,800))    
V = [Point{3,Float64}(x[i], y[i], z[i]) for i in eachindex(x)]

function toquadmesh(V, numRows)
    numCols = length(V)/numRows # Number of point "columns"
        if !isinteger(numCols) || numCols<1
        throw(ArgumentError("The length(V)/numRows should produce an integer >1 but is instead $numCols"))
    end
    numCols = Int(numCols)
    ij2ind(i,j) = i + numCols*(j-1) # function to convert subscript to linear indices
    F = Vector{QuadFace{Int}}(undef,(numCols-1)*(numRows-1))     
    i_f = 1     
    @inbounds for i = 1:(numCols-1)
        @inbounds for j = 1:(numRows-1)
            F[i_f] = QuadFace{Int}(ij2ind(i,j),
                                   ij2ind(mod1(i+1,numCols), j),
                                   ij2ind(mod1(i+1,numCols), mod1(j+1,numRows)), 
                                   ij2ind(i,mod1(j+1,numRows)))
            i_f += 1            
        end
    end
    return F
end

for testCase = 1:1:length(methodSet)
    methodNow = methodSet[testCase]
    z_i = itp(x_i,y_i; method=methodNow)
    V_i = [Point{3,Float64}(x_i[i], y_i[i], z_i[i]) for i in eachindex(x_i)]   
    if testCase<=4
        r = 1
        c = testCase
    else
        r = 2
        c = testCase-4
    end
    F_i = toquadmesh(V_i, numPointsInterp)
    
    ax = Axis3(fig[r, c], clip=false, viewmode = :free, aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title="method: " * titleStrings[testCase])
    hp1 = scatter!(ax, V, markersize=12, color=:black)
    # hp2 = scatter!(ax, V_i, markersize=8, color=:red)
    hp3 = mesh!(ax, GLMakie.GeometryBasics.Mesh(V_i,F_i), color=z_i)
end
fig

[DanielVandH]

Yeah, unfortunately extrapolation outside of the convex hull is a complex issue. How extrapolation is implemented is here https://github.com/DanielVandH/NaturalNeighbours.jl/blob/main/src/interpolation/extrapolation.jl, where I am basically just projecting the point onto the convex hull and then using the data on that edge to interpolate.

I did originally want to implement better extrapolation techniques - see the fourth point here #1, it's called "ghost point extrapolation". That would be something good to aim for based on the result in Bobach's thesis as well, but I was never able to really understand their papers and (frustratingly...) they provide no code. Maybe there is something simpler to be done in the interim, but I don't have a lot of time for this package these days to do it myself sadly.

[Kevin-Mattheus-Moerman]

@DanielVandH thanks for the quick reply Daniel.

I understand extrapolation outside the convex hull is tricky. And also of course I understand your time constraints.

I'll just add that it does look like the non-smoothness starts inside the convex hull. In the video below I'm only interpolating inside the convex hull only, but you can see that there are "origami folds" at the boundary. Can this be avoided? Perhaps this is because the input is sparse at the boundary i.e. not as well sampled as the interior yet the function does have a big difference from one point to the next.

naturalNeighbourInterp_02.mp4

If I give it a domain where interior and boundary are well sampled it looks more well behaved.

batman_interp.mp4

[DanielVandH]

Can you also show a plot of the Voronoi tessellation with the convex hull plotted on top?