Interpolations extension

This extension overloads the interpolate function from Interpolations.jl to provide a way to construct a global interpolant from the discrete data in a LevelSet or LevelSetEquation. This can be useful in situations where you want to evaluate the approximate underlying functions at points that are not on the grid.

Here is an example of how to construct such an interpolant:

using LevelSetMethods, Interpolations
LevelSetMethods.set_makie_theme!()
a, b = (-2, -2), (2, 2)
ϕ   = LevelSetMethods.star(CartesianGrid(a, b, (50, 50)))
itp = interpolate(ϕ, BSpline(Cubic())) # create the interpolant
50×50 scale(interpolate(OffsetArray(::Matrix{Float64}, 0:51, 0:51), BSpline(Cubic(Line(OnGrid())))), (LinRange{Float64}(-2.0, 2.0, 50), LinRange{Float64}(-2.0, 2.0, 50))) with element type Float64:
 1.65165  1.57711  1.50533  1.43686  …  1.43686  1.50533  1.57711  1.65165
 1.61386  1.5362   1.46098  1.38873     1.38873  1.46098  1.5362   1.61386
 1.57995  1.49933  1.42076  1.34479     1.34479  1.42076  1.49933  1.57995
 1.54996  1.4666   1.38487  1.30531     1.30531  1.38487  1.4666   1.54996
 1.52391  1.43809  1.35346  1.27052     1.27052  1.35346  1.43809  1.52391
 1.50169  1.41378  1.32658  1.24056  …  1.24056  1.32658  1.41378  1.50169
 1.48312  1.39355  1.30421  1.2155      1.2155   1.30421  1.39355  1.48312
 1.46789  1.37715  1.28617  1.19527     1.19527  1.28617  1.37715  1.46789
 1.45557  1.36419  1.27215  1.17966     1.17966  1.27215  1.36419  1.45557
 1.44556  1.35414  1.26166  1.16825     1.16825  1.26166  1.35414  1.44556
 ⋮                                   ⋱                             
 1.36698  1.32382  1.28322  1.2451      1.2451   1.28322  1.32382  1.36698
 1.44772  1.40655  1.36761  1.33074     1.33074  1.36761  1.40655  1.44772
 1.52915  1.48942  1.45155  1.41531     1.41531  1.45155  1.48942  1.52915
 1.61063  1.57183  1.5345   1.49835     1.49835  1.5345   1.57183  1.61063
 1.69165  1.65333  1.61608  1.57959  …  1.57959  1.61608  1.65333  1.69165
 1.77183  1.7336   1.69606  1.65887     1.65887  1.69606  1.7336   1.77183
 1.85092  1.81245  1.77431  1.73615     1.73615  1.77431  1.81245  1.85092
 1.92875  1.88976  1.85079  1.81147     1.81147  1.85079  1.88976  1.92875
 2.0052   1.9655   1.92554  1.88494     1.88494  1.92554  1.9655   2.0052

Once constructed, the interpolant can be used to evaluate the level-set function anywhere inside the grid:

itp(0.5, 0.5)
-0.11611218885936599

This can be used e.g. to plot the level-set function using Makie:

using GLMakie
xx = yy = -2:0.01:2
contour(xx, yy, [itp(x,y) for x in xx, y in yy]; levels = [0], linewidth = 2)
Example block output

Trying to evaluate it outside the domain will throw an error:

try
  itp(3, 0.1)
catch e
    println("Error caught")
end
Error caught

Using it on three-dimensional level sets is similar:

using LinearAlgebra
grid = CartesianGrid((-1.5, -1.5, -1.5), (1.5, 1.5, 1.5), (50, 50, 50))
P1, P2 = (-1, 0, 0), (1, 0, 0)
b = 1.05
f = (x) -> norm(x .- P1)*norm(x .- P2) - b^2
ϕ = LevelSet(f, grid)
itp = interpolate(ϕ) # cubic spline by default
println("ϕ(0.5, 0.5, 0.5)   = ", f((0.5, 0.5, 0.5)))
println("itp(0.5, 0.5, 0.5) = ", itp(0.5, 0.5, 0.5))
ϕ(0.5, 0.5, 0.5)   = 0.3336406616345071
itp(0.5, 0.5, 0.5) = 0.33364017975248