Interpolation

LevelSetMethods.jl provides a built-in high-performance piecewise polynomial interpolation scheme. This allows you to construct a continuous interpolant from the discrete data in a LevelSet or LevelSetEquation. This is useful for evaluating the level-set function at arbitrary coordinates, computing analytical gradients and Hessians, or for visualization.

Basic Usage

To construct an interpolant, use the interpolate function:

using LevelSetMethods
a, b = (-2.0, -2.0), (2.0, 2.0)
ϕ   = LevelSetMethods.star(CartesianGrid(a, b, (50, 50)))
# Add boundary conditions for safe evaluation near edges
bc = ntuple(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), 2)
ϕ = LevelSetMethods.add_boundary_conditions(ϕ, bc)

itp = interpolate(ϕ) # cubic interpolation by default (order=3)
PiecewisePolynomialInterpolant
  ├─ order: 3
  └─ field: MeshField on CartesianGrid in ℝ²
     ├─ domain:  [-2.0, 2.0] × [-2.0, 2.0]
     ├─ nodes:   50 × 50
     ├─ spacing: h = (0.08163, 0.08163)
     ├─ bc:     Degree 1 extrapolation (all)
     ├─ eltype:  Float64
     └─ values:  min = -1.121,  max = 2.005

The returned object is a PiecewisePolynomialInterpolant, which is callable and efficient. Once constructed, the interpolant can be used to evaluate the level-set function anywhere inside (and even slightly outside, using boundary conditions) the grid:

itp(0.5, 0.5)
-0.11616315120712294

Plotting

Interpolation is particularly useful for creating smooth plots. Here is an example using Makie:

using GLMakie
LevelSetMethods.set_makie_theme!()
xx = yy = -2:0.01:2
# Evaluate on a fine grid for plotting
contour(xx, yy, [itp(x, y) for x in xx, y in yy]; levels = [0], linewidth = 2)
Example block output

Derivatives

The interpolant also provides analytical gradients and Hessians via the LevelSetMethods.gradient and LevelSetMethods.hessian functions. These are zero-allocation and computed using Horner's method on the underlying Lagrange polynomials.

x = (0.1, 0.2)
val  = itp(x)
grad = LevelSetMethods.gradient(itp, x)
hess = LevelSetMethods.hessian(itp, x)
println("Value:    ", val)
println("Gradient: ", grad)
Value:    -0.9335268915671028
Gradient: [3.0106203502669895, -1.5126156350576423]

Three-dimensional Level Sets

Using it on three-dimensional level sets is identical:

using LinearAlgebra, StaticArrays
grid = CartesianGrid((-1.5, -1.5, -1.5), (1.5, 1.5, 1.5), (32, 32, 32))
P1, P2 = (-1.0, 0.0, 0.0), (1.0, 0.0, 0.0)
b = 1.05
f = (x) -> norm(x .- P1)*norm(x .- P2) - b^2
ϕ3 = LevelSet(f, grid)
bc3 = ntuple(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), 3)
ϕ3 = LevelSetMethods.add_boundary_conditions(ϕ3, bc3)

itp3 = interpolate(ϕ3)
println("ϕ(0.5, 0.5, 0.5)   = ", f(SVector(0.5, 0.5, 0.5)))
println("itp(0.5, 0.5, 0.5) = ", itp3(0.5, 0.5, 0.5))
ϕ(0.5, 0.5, 0.5)   = 0.3336406616345071
itp(0.5, 0.5, 0.5) = 0.33361166021372846