Spline Interpolation Model: Methods & Applications

• Spline interpolation models are defined by constructing piecewise low-degree functions with continuity constraints on derivatives at the knots.
• They employ computational strategies such as adaptive knot allocation and tridiagonal solves to optimize both local accuracy and overall efficiency.
• These models extend from classical cubic splines to variational, geometric, and trainable frameworks, finding applications in signal processing, neuroimaging, and manifold interpolation.

Spline interpolation model denotes a class of interpolation constructions in which prescribed data are matched by piecewise low-degree functions, or by generalized basis functions that play the same role, joined at knots under explicit smoothness constraints. In the material considered here, the model ranges from classical piecewise-cubic forms and B-spline bases to $q$-polynomial, rational-fractal, manifold, Wasserstein, and image-space generalizations; across these variants, the recurring objective is to satisfy interpolation conditions while controlling curvature, acceleration, or an analogous roughness quantity adapted to the ambient domain (Herscovici, 2018, Sire et al., 13 Oct 2025, Justiniano et al., 2023).

1. Core definition and structural invariants

A representative classical specification writes a cubic spline as a piecewise cubic function

$$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$

with interpolation constraints

$S(x_i)=y_i,\qquad i=0,1,\ldots,n,$

and continuity at internal knots of the function, first derivative, and second derivative: $$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$ One formulation in the literature further imposes periodic endpoint matching,

$$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$

while other models use clamped, not-a-knot, or Hermite endpoint conditions (Yu, 2023).

The invariants that persist across the surveyed models are more fundamental than any single boundary prescription. First, spline interpolation is local in the sense that the interpolant is assembled interval by interval, or from compactly supported basis functions. Second, smoothness is encoded by continuity of derivatives or by an intrinsic analogue such as covariant derivatives, Laplace–Beltrami curvature, or Wasserstein acceleration. Third, the interpolation conditions may act on function values alone, or on Hermite data that include derivatives at the knots. The bounded-interval local scheme, for example, requires exact interpolation at the sample sites together with endpoint Hermite matching up to order $m-1$ and polynomial reproduction up to degree $m-1$ (Walt, 2015).

A central source of variation is the endpoint model. The cubic $q$-spline replaces ordinary endpoint slopes by clamped Jackson $q$-derivatives,

$$D_q S_1(x_0)=D_q f(x_0),\qquad D_q S_n(x_n)=D_q f(x_n),$$

whereas AutoKnots uses a cubic spline with not-a-knot endpoint conditions,

$$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$0

and the real-time local bounded-interval construction enforces endpoint derivative matching through dedicated local molecules (Herscovici, 2018, Vitenti et al., 2024, Walt, 2015). Accordingly, “spline interpolation model” is best understood not as a single formula but as a structural template whose concrete realization is fixed by basis choice, knot placement, smoothness operator, and boundary conditions.

2. Knots, support, and computational structure

Knot placement governs both approximation quality and computational cost. In the bounded-interval local construction, the knot sequence is chosen from sample midpoints,

$$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$1

and the final interpolant is formed by a blending operator

$$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$2

Here $$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$3 supplies local quasi-interpolation and polynomial reproduction, whereas $$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$4 corrects the residual so that exact interpolation and endpoint Hermite conditions hold. The significance of this formulation is that it avoids the global linear solves associated with classical spline interpolation while preserving compact support and optimal interior approximation order $$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$5 (Walt, 2015).

Adaptive knot allocation changes the knot set itself from a fixed design variable into an error-controlled output. AutoKnots begins with 6 uniformly spaced knots on $$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$6, constructs a not-a-knot cubic spline, and then repeatedly inserts interval midpoints into the knot set

$$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$7

whenever local tests are not yet certified. Acceptance is based on a pointwise midpoint criterion and an integral criterion using Simpson’s $$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$8 rule: $$S(x)= \left\{ S_i(x)=a_i x^3+b_i x^2+c_i x+d_i,\quad x\in [x_i,x_{i+1}],\ i=0,1,\ldots,n-1 \right\},$$9 Every newly evaluated point becomes a permanent knot, so the adaptive search directly constructs the final interpolation table. The reported motivation is minimal manual configuration with a user-specified precision target, together with nonuniform allocation of knots to difficult regions (Vitenti et al., 2024).

A different computational motif is the tridiagonal solve. In the cubic $S(x_i)=y_i,\qquad i=0,1,\ldots,n,$0-spline, the unknown nodal moments $S(x_i)=y_i,\qquad i=0,1,\ldots,n,$1 satisfy a linear system

$S(x_i)=y_i,\qquad i=0,1,\ldots,n,$2

with tridiagonal matrix structure. This preserves the familiar cubic-spline computational pattern while replacing ordinary polynomials by $S(x_i)=y_i,\qquad i=0,1,\ldots,n,$3-polynomials and ordinary derivatives by Jackson $S(x_i)=y_i,\qquad i=0,1,\ldots,n,$4-derivatives (Herscovici, 2018). The same preference for low-bandwidth linear algebra reappears in the AutoKnots implementation, where the not-a-knot cubic spline is reduced to a symmetric positive-definite tridiagonal system solved by LAPACK dptsv (Vitenti et al., 2024).

These constructions show that spline interpolation models differ not only in analytic definition but in their computational semantics. Some are global variational solvers; some are local evaluation operators; some are adaptive knot generators. This suggests that “model” includes the interpolation space and smoothness constraints, but also the algorithmic pathway by which the interpolant is made numerically available.

3. Variational and geometric spline models

In non-Euclidean settings, the spline model is no longer most naturally defined by a piecewise polynomial ansatz. Instead, the model is specified by an intrinsic energy whose minimizers play the role of splines (Sire et al., 13 Oct 2025, Lin et al., 2023, Justiniano et al., 2021, Justiniano et al., 2023, Chewi et al., 2020).

Model family	Defining smoothness quantity	Domain
Manifold thin-plate spline	$S(x_i)=y_i,\qquad i=0,1,\ldots,n,$5	Compact Riemannian manifold
Higher-order Riemannian spline	$S(x_i)=y_i,\qquad i=0,1,\ldots,n,$6	Riemannian manifolds and Lie groups
Metamorphosis spline	Flow acceleration + second material derivative	Space of images
Wasserstein variational spline	Wasserstein acceleration regularized by action	Probability measures
Transport spline	OT coupling + Euclidean natural cubic spline	$S(x_i)=y_i,\qquad i=0,1,\ldots,n,$7

On a compact Riemannian manifold $S(x_i)=y_i,\qquad i=0,1,\ldots,n,$8, the spline interpolant is the minimizer of

$S(x_i)=y_i,\qquad i=0,1,\ldots,n,$9

subject to $$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$0, or, in the smoothing case,

$$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$1

The model is expressed in RKHS form through kernels built from the Laplace–Beltrami spectrum, but the paper’s main contribution is computational: spline interpolation is identified with kriging for a Gaussian field, then approximated by finite elements on a triangulated mesh and implemented via sparse intrinsic GMRF precision matrices. The same framework admits anisotropy by deforming the Riemannian metric, thereby changing the Laplace–Beltrami operator itself rather than merely postulating an anisotropic covariance (Sire et al., 13 Oct 2025).

Higher-order Riemannian spline interpolation generalizes the Euclidean squared-acceleration principle by replacing ordinary derivatives with covariant derivatives along a curve. The energy

$$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$2

leads to the Euler–Lagrange operator

$$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$3

together with higher-order concurrency conditions at interior knots. The paper studies both exact interpolation and least-squares fitting through successively connected curve networks, and solves the problem by gradient flows whose asymptotic stationary limits furnish spline solutions (Lin et al., 2023).

For image-valued data, the metamorphosis model introduces a spline functional that combines the Eulerian acceleration of the flow and the second material derivative of image intensity. In its regularized relaxed form,

$$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$4

This is an image-space analogue of cubic spline interpolation rather than a piecewise geodesic interpolant. A variational time discretization is shown to converge to the relaxed continuous model via Mosco convergence, and the fully discrete implementation uses finite differences together with stable cubic B-spline interpolation of deformed quantities (Justiniano et al., 2021).

Wasserstein-space splines appear in two distinct forms. A variational formulation defines a continuous spline energy

$$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$5

and regularizes it by the usual action functional. Its time-discrete counterpart replaces Euclidean second differences by Wasserstein barycentric deviations and is computed numerically with entropy-regularized OT and Sinkhorn iterations (Justiniano et al., 2023). A different construction, the transport spline, first couples measures by successive Monge maps and then fits a natural cubic spline in Euclidean space to the induced particle trajectories. In one dimension, when the quantile spline remains a valid quantile function, this transport spline coincides with the E-spline (Chewi et al., 2020).

4. Generalized bases, shape control, and fractalization

Several models generalize the polynomial basis itself rather than only the ambient geometry. The cubic $$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$6-spline is the clearest algebraic example. It is a $$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$7-analogue of the classical cubic spline based on Jackson’s $$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$8-derivative

$$S_i(x_{i+1})=S_{i+1}(x_{i+1}),\qquad S_i^{(1)}(x_{i+1})=S_{i+1}^{(1)}(x_{i+1}),\qquad S_i^{(2)}(x_{i+1})=S_{i+1}^{(2)}(x_{i+1}).$$9

with piecewise segments that are $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$0-polynomials of degree at most three. The construction enforces continuity of $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$1, $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$2, and $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$3, and each segment is expressed explicitly in terms of $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$4-polynomials $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$5. As $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$6, one has

$$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$7

so the model tends to the standard cubic spline. The parameter $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$8 is therefore not purely formal; in the reported example for $$S^{(1)}(a+0)=S^{(1)}(b-0),\qquad S^{(2)}(a+0)=S^{(2)}(b-0),$$9, decreasing $m-1$0 intensifies oscillation while increasing $m-1$1 reduces it and gives “significant improvement” (Herscovici, 2018).

Domain-informed spline interpolation replaces shift-invariant B-spline bases by shift-variant generators adapted to a known inhomogeneous domain. The model assumes nonnegative subdomain descriptor functions $m-1$2 satisfying

$m-1$3

and defines local basis functions $m-1$4 by decomposing each basis into dominant and residual parts. The resulting interpolation space

$m-1$5

preserves partition of unity and the Riesz basis property, and reduces to the ordinary shifted B-spline basis on homogeneous regions. The practical implication is that neighboring samples are weighted not only by distance but by domain similarity, which was demonstrated in Monte Carlo simulations and neuroimaging interpolation guided by tissue-probability maps (Behjat et al., 2018).

Fractalization introduces self-reference into the spline model. The $m-1$6 shape-preserving rational cubic spline FIF is defined by the functional equation

$m-1$7

where $m-1$8 are scaling factors and $m-1$9 are tension or shape parameters appearing in

$m-1$0

When $m-1$1, the model reduces to the classical Delbourgo–Gregory rational cubic spline; with nonzero $m-1$2, it becomes a fractal generalization with explicit sufficient conditions for positivity, monotonicity, and convexity (Chand et al., 2015).

The cubic spline super fractal interpolation function extends the same idea from a single IFS to a super IFS, allowing the graph to be generated by different local maps at different scales. Under equidistant partitions and suitable smallness conditions on the vertical scaling parameters, the interpolant and its first two derivatives satisfy

$m-1$3

This preserves spline-like convergence order while permitting multiscale variability that a single classical cubic spline does not model (Kapoor et al., 2012).

5. Optimization, kernel design, and trainable spline models

A distinct line of work treats the spline basis itself as an optimization variable. In optimized spline interpolation, one searches over a compact-support class

$m-1$4

for a basis $m-1$5 whose induced interpolating kernel

$m-1$6

best approximates either a target signal $m-1$7 or a desired interpolation filter $m-1$8 in $m-1$9. The variational normal equation is reduced to a finite Hermitian linear system on the support interval, so the design problem becomes tractable while preserving compact support. For sinc approximation, the reported SNR values are 20.39 dB for the optimized spline and 13.15 dB for the cubic spline, and image-resizing experiments report average PSNR 28.88 dB for the lowpass-optimized spline in one setting (Madani et al., 2011, Madani et al., 2010).

Smooth supersaturated models arrive at a related objective from an algebraic direction. One deliberately extends an interpolating monomial basis beyond the saturated case and then minimizes a roughness functional such as

$q$0

over the free coefficients while maintaining exact interpolation. In one dimension the resulting interpolants approach cubic-spline smoothness rapidly as extra basis terms are added; in two dimensions they act as an alternative to thin-plate splines (Bates et al., 2008). This suggests that spline behavior can emerge from a smoothness-constrained basis extension even when the basis is global and polynomial rather than piecewise polynomial.

The zero-delay smoothing spline model reinterprets online spline construction as sequential decision making. After each incoming sample $q$1, the method must irrevocably choose a polynomial segment

$q$2

subject to continuity constraints encoded by inherited boundary derivatives. The per-step cost is

$q$3

and the proposed policy augments this myopic objective by a learned cost-to-go term generated by a recurrent neural network. The resulting architecture is a recurrent network plus a differentiable convex optimization layer, so continuity constraints remain exact while the policy becomes future-aware (Ruiz-Moreno et al., 2022).

Taken together, these models shift spline interpolation from a fixed analytic family to an optimization framework. This suggests that, in contemporary usage, a spline interpolation model may be selected not only by polynomial order or knot sequence but by the target filter, data source, roughness penalty, or causal deployment constraint.

6. Application patterns, evaluation, and recurrent limitations

The application literature shows that spline interpolation often functions as an internal module rather than an end in itself. In stock-market forecasting, cubic spline interpolation is explicitly used as a preprocessing layer that fills missing values before stationarity testing, differencing, and ARIMA fitting. The reported pipeline is: original series with gaps, cubic spline interpolation, ADF testing, differencing, ARIMA order selection, fitting, and forecasting. The paper states that spline interpolation is not presented as a direct forecasting model, and it does not provide an ablation against raw ARIMA without interpolation. Its rolling-window results show MSE around 10 when test size is 10, rising to around 700 when test size is 1000, supporting the claim that the method is suitable mainly for short-term forecasting and is compromised by intense fluctuations (Yu, 2023).

In CNC trajectory generation, the spline interpolation model is used for error-controlled $q$4-continuous tool paths. The proposed cubic B-spline uses, on each interval, endpoint matching of the $q$5-, $q$6-, and $q$7-order derivatives, achieved by introducing two internal knots rather than increasing polynomial degree. The local error satisfies

$q$8

which yields an adaptive step-size rule

$q$9

for prescribed tolerance $q$0. The paper also notes that for straight-line segments the corresponding $q$1 is 0, so no maximal step-length restriction is imposed there (Ze-Wei et al., 2024).

In animation line inbetweening, thin-plate splines provide a coarse nonrigid motion estimate from matched keypoints between two raster line drawings. The TPS transformation uses the standard kernel

$q$2

and bending energy

$q$3

A lightweight motion-refinement module and a UNet then convert the coarse warp into the final intermediate frame. The reported user study assigns 70.50% overall preference to the proposed method, compared with 19.40% and 10.10% for two baselines, and the proposed Weighted Chamfer Distance achieves 72.50% consistency with human preference across settings (Zhu et al., 2024).

Other application domains broaden the notion of what is being interpolated. Domain-informed B-spline interpolation is used in neuroimaging, where anatomical segmentations guide interpolation of functional maps (Behjat et al., 2018). AutoKnots is integrated into NumCosmo and used to interpolate costly halo-profile quantities such as $q$4 and $q$5 while meeting prescribed tolerances (Vitenti et al., 2024). Manifold spline interpolation is validated on the sphere through analytical test cases and a pollution-related study, and also on the surface of a cylinder (Sire et al., 13 Oct 2025). Wasserstein splines are applied to synthesized textures by interpolating feature distributions rather than pixels (Justiniano et al., 2023).

Several recurrent limitations are also explicit in the literature. Local midpoint-based adaptivity can miss plateau-like or low-curvature regions unless a secondary refinement heuristic is used (Vitenti et al., 2024). Wasserstein spline interpolants may exhibit overshooting, in close analogy with classical cubic-spline overshoot (Justiniano et al., 2023). Hybrid forecasting pipelines may rely on spline interpolation without isolating its quantitative gain through a direct baseline comparison (Yu, 2023). More generally, the surveyed papers show that spline interpolation is not synonymous with a universal cubic recipe: it may be exact or smoothing, local or global, deterministic or learned, Euclidean or geometric, and the adequacy of a given model depends on whether the underlying notion of smoothness is polynomial, variational, transport-based, or shape constrained.