Shape-preserving interpolation on simplicial meshes

When to interpolate instead of fitting

fit estimates one global polynomial (or one per glued cell) by least squares. When the response is piecewise smooth, with kinks where an active set or a regime changes, a global polynomial of moderate degree oscillates between design points and can overshoot the data. If the response is optimized afterwards, an overshoot becomes a spurious optimum that the data never showed.

hermite_interpolant() builds a local alternative from data at the vertices of a simplicial mesh: a continuous piecewise cubic Bézier model that passes through every datum and, with the default clamp=True, never leaves the range of the data at the vertices of the simplex it is evaluated in. On a chain of intervals it is the familiar monotone piecewise cubic Hermite interpolation (PCHIP) [FC80]; on triangles and tetrahedra it is its natural generalization, which we call simplicial PCHIP.

Use it when evaluations are exact (or nearly so) and the design points are given, for example a grid or the points chosen by an adaptive loop. With noisy observations prefer a regularized fit: interpolation reproduces the noise.

Interpolate on a mesh

A SimplicialMesh takes vertex coordinates of shape (n_vertices, d) and simplices of shape (n_simplices, d + 1) given by vertex numbers. The simplices must not be degenerate and must meet face to face. The example triangulates a jittered grid of the unit square and interpolates a function with a kink along a diagonal.

import torch
from torch_bsf.interpolation import SimplicialMesh, hermite_interpolant

n = 6
i, j = torch.meshgrid(torch.arange(n + 1), torch.arange(n + 1), indexing="ij")
vertices = torch.stack([i.flatten(), j.flatten()], dim=1).double() / n
number = lambda a, b: a * (n + 1) + b
simplices = []
for a in range(n):
    for b in range(n):
        simplices.append([number(a, b), number(a + 1, b), number(a + 1, b + 1)])
        simplices.append([number(a, b), number(a, b + 1), number(a + 1, b + 1)])
mesh = SimplicialMesh(vertices, simplices)

def response(x):
    return torch.stack([(x[:, 0] - x[:, 1]).abs() + x[:, 0] ** 2, x.sum(1)], dim=1)

model = hermite_interpolant(mesh, response(vertices))
x = torch.rand(1000, 2, dtype=torch.float64)
points = mesh.locate(x)
prediction = model(points)
assert torch.allclose(model(mesh.locate(vertices)), response(vertices), atol=1e-12)
assert prediction.shape == (1000, 2)

mesh.locate returns Points: the cell ID of the simplex containing each point and its barycentric coordinates there. model is an ordinary BezierModel on mesh.space, a BezierSpace of degree three compiled from the mesh. Everything that accepts a model, such as differentiation with respect to the points or subdivide(), accepts this one.

The barycentric coordinates of cell mesh.cell_ids[k] refer to the vertices mesh.simplices[k], which lists the vertex numbers of simplex k in ascending order (the order you passed is not kept). A point on a face shared by several simplices is assigned to one of them; the prediction does not depend on the choice because the interpolant is continuous. Points outside the mesh by more than tol raise ValueError, while points within tol of its boundary are projected onto it. Point location tests every simplex, which is simple and exact but costs O(n_points * n_simplices).

How the cubic pieces are built

On each simplex the model is a cubic Bézier simplex. Its controls are set from the data \(f_i\) and gradients \(g_i\) at the vertices \(x_i\), like the Bézier form of a cubic Hermite segment:

  • at vertex \(i\), the control is \(f_i\);

  • on edge \(i \to j\), next to vertex \(i\), it is \(f_i + g_i \cdot (x_j - x_i) / 3\);

  • at the centre of the triangle face \((i, j, l)\), it is the sum of the six edge controls of that face divided by four, minus \((f_i + f_j + f_l) / 6\).

A cubic simplex has no other controls. The face formula is the degree elevation of the quadratic Bézier triangle whose edge midpoint controls match the cubic edge controls [Far02], so the interpolant reproduces every quadratic polynomial when the gradients are exact.

Each control depends only on the data of the edge or face that carries it. SimplicialMesh glues simplices that share a face so that the face’s controls are the same coefficient rows from both sides, which makes the interpolant continuous (C0). It is not C1 in general: gradients across a shared face can jump.

When gradients is omitted, estimate_gradients() fits a quadratic by weighted least squares (weights \(1/|x_w - x_v|^2\)) to the data on the smallest \(k\)-ring of mesh neighbours, \(k \ge 2\), that has two points more than the quadratic has unknowns, and falls back to a linear fit on the one-ring. Pass your own gradients of shape (n_vertices, d, n_outputs) if you know them.

from torch_bsf.interpolation import estimate_gradients

quadratic = 1 + vertices[:, :1] - vertices[:, 1:] ** 2 + vertices[:, :1] * vertices[:, 1:]
exact = hermite_interpolant(mesh, quadratic, clamp=False)
truth = 1 + x[:, :1] - x[:, 1:] ** 2 + x[:, :1] * x[:, 1:]
assert torch.allclose(exact(points), truth, atol=1e-10)
assert estimate_gradients(mesh, quadratic).shape == (mesh.n_vertices, 2, 1)

Range limiting and its error guarantee

With clamp=True (the default) each edge control is clipped, output by output, into the range of the data at the two vertices of its edge, and each face control into the range of the data at the three vertices of its face. Clipping uses only the data of that edge or face, so shared controls are still shared and continuity is kept. Vertex data are unchanged.

Bernstein polynomials are nonnegative and sum to one, so a Bézier simplex takes values in the convex hull of its controls. Every control of a simplex \(T\) is now within the range of its vertex data, hence for every output

\[\min_{v \in T} f(v) \le p(x) \le \max_{v \in T} f(v), \qquad x \in T.\]

If the data come from a function \(f\), then for \(x \in T\),

\[|p(x) - f(x)| \le \max_{v \in T} |f(v) - f(x)| \le \omega_f(h),\]

where \(h\) is the longest edge of the mesh and \(\omega_f(\delta) = \sup\{|f(y) - f(z)| : |y - z| \le \delta\}\) is the modulus of continuity of \(f\). The bound holds for any gradients, including poorly estimated ones, and for kinked or discontinuous-derivative responses; refining the mesh therefore converges uniformly for any continuous \(f\) on the meshed domain. In particular, the interpolant never predicts a value outside the observed range of a cell, so its minimum over the domain equals the smallest datum and is attained at a design point: minimizing the interpolant cannot produce a spurious optimum.

y = response(vertices)
corners = y[mesh.simplices]
cell_of = {c: k for k, c in enumerate(mesh.cell_ids.tolist())}
rows = torch.tensor([cell_of[c] for c in points.cell_ids.tolist()])
assert (prediction >= corners.amin(dim=1)[rows] - 1e-12).all()
assert (prediction <= corners.amax(dim=1)[rows] + 1e-12).all()

In one dimension, clipping keeps the slope ratios \(\alpha = g_i h / \Delta\) and \(\beta = g_{i+1} h / \Delta\) of each interval in \([0, 3]\), where \(\Delta\) is the change of the data over the interval. This is the box condition of Fritsch and Carlson, which makes the cubic monotone on every interval where the data are monotone. Without clipping, a step in the data overshoots:

line = torch.arange(8, dtype=torch.float64)[:, None]
chain = SimplicialMesh(line, [[k, k + 1] for k in range(7)])
step = torch.tensor([0., 0., 0., 0., 1., 1., 1., 1.], dtype=torch.float64)
t = chain.locate(torch.linspace(0, 7, 701, dtype=torch.float64)[:, None])
loose = hermite_interpolant(chain, step, clamp=False)(t)
tight = hermite_interpolant(chain, step)(t)
print(f"unclamped: [{loose.min():.3f}, {loose.max():.3f}]")
print(f"clamped:   [{tight.min():.3f}, {tight.max():.3f}]")
assert (tight.diff(dim=0) >= -1e-12).all()

unclamped: [-0.047, 1.047]
clamped:   [0.000, 1.000]

Clipping trades accuracy for safety. Near a local extremum of the response, the true function leaves the range of the vertex data of a cell, and the clamped interpolant cannot follow it, so quadratic data are no longer reproduced exactly in such cells. Use clamp=False when the response is smooth and exact quadratic reproduction matters more than the range guarantee.

Choose where to add design points

refinement_indicator() scores each simplex by its volume times the difference between the cubic interpolant and linear interpolation at its centroid. A large score means the data bend within a large cell. Evaluating the response near the midpoint of the longest edge of the top-scoring cell and remeshing is a simple adaptive design. The example refines a one-dimensional kink, where remeshing is just sorting.

from torch_bsf.interpolation import refinement_indicator

def kink(s):
    return (s - 0.37).abs()

def chain_mesh(s):
    s = s.sort().values
    return s, SimplicialMesh(s[:, None], [[k, k + 1] for k in range(len(s) - 1)])

test = torch.linspace(0, 1, 1001, dtype=torch.float64)
design, chain = chain_mesh(torch.linspace(0, 1, 5, dtype=torch.float64))
first = (hermite_interpolant(chain, kink(design))(chain.locate(test[:, None]))[:, 0] - kink(test)).abs().max()
for _ in range(6):
    worst = int(refinement_indicator(chain, kink(design)).argmax())
    a, b = chain.simplices[worst].tolist()
    design, chain = chain_mesh(torch.cat([design, (design[a:a + 1] + design[b:b + 1]) / 2]))
last = (hermite_interpolant(chain, kink(design))(chain.locate(test[:, None]))[:, 0] - kink(test)).abs().max()
assert last < first / 4

In higher dimensions, LongestEdge in torch_bsf.splitting selects the longest edge of a cell, and a triangulation library rebuilds the mesh. With the optional SciPy dependency (the sparse extra), a Delaunay triangulation can be passed directly:

from scipy.spatial import Delaunay

scattered = torch.cat([torch.tensor([[0., 0.], [1., 0.], [0., 1.], [1., 1.]]), torch.rand(30, 2)]).double()
delaunay = SimplicialMesh(scattered, Delaunay(scattered.numpy()).simplices)
smooth = hermite_interpolant(delaunay, response(scattered))
score = refinement_indicator(delaunay, response(scattered))
assert score.shape == (delaunay.n_simplices,)

Limitations

  • The mesh is not generated here; the simplices must be full-dimensional, non-degenerate and conforming (two simplices meet in a common face or not at all). Otherwise shared controls are not shared and the interpolant is discontinuous.

  • The interpolant is cubic and C0. C1 macro-elements (Clough–Tocher, Powell–Sabin) and non-simplicial cells such as prisms are not provided.

  • Quadratic reproduction requires enough neighbours around each vertex for the gradient estimate; near a sparse boundary only linear data are reproduced unless you pass gradients.

  • Computation is in float64 on the CPU.