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)

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