Interpolation — Polynomial, Spline & B-spline Interpolation

Overview

The sangi interpolation module provides a family of functions that reconstruct intermediate values from known sample points $(x_i, y_i)$. It is organized into the following categories by use case.

• Piecewise & polynomial interpolation — Lagrange / Newton divided differences / linear & bilinear / zero-order hold (step) / barycentric Lagrange
• Cubic spline — $C^2$-continuous piecewise cubic curves with natural or clamped boundaries
• Shape-preserving & local interpolation — PCHIP (monotonicity-preserving), modified Akima, Catmull-Rom, smooth interpolation, Floater-Hormann rational interpolation, quintic Hermite
• B-spline — Cox-de Boor basis, clamped knots, de Boor evaluation, B-spline interpolation / regression
• Trigonometric interpolation — trigonometric-polynomial interpolation of equally spaced periodic data
• Spline fitting — least-squares fits with a smoothing parameter (1D / parametric / 2D tensor product)
• Image interpolation — nearest-neighbor, bicubic (Catmull-Rom kernel), and Lanczos resampling

The template parameter $T$ requires concepts::Field or concepts::OrderedField. Typically you pass double; when higher precision is needed you can pass the arbitrary-precision type Float. Everything lives in namespace sangi.

Many functions provide both an overload taking std::span<const T> and one taking std::vector<T>. This page mainly describes the std::vector overloads (the arguments of the span overloads have the same meaning).

Build

// Piecewise polynomial, spline, shape-preserving, fitting, image interpolation
#include <math/interpolation/Interpolation.hpp>

// De Boor basis, cubicSpline, B-spline regression
#include <math/interpolation/BSpline.hpp>

Header-only. No library linkage required. The functions that solve a linear system internally (B-spline interpolation, spline fitting, surface spline) depend on <math/linalg/solvers.hpp>, but Interpolation.hpp includes it automatically.

Piecewise & Polynomial Interpolation

Function reference

lagrange_interpolation

template<concepts::Field T>
T lagrange_interpolation(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T xi);

Behaviour: directly evaluates, at the point $x_i$, the degree-$(n-1)$ Lagrange interpolating polynomial $L(x) = \sum_i y_i \prod_{j \ne i} \frac{x - x_j}{x_i - x_j}$ through $n$ sample points. Each evaluation is $O(n^2)$.

Parameters:

Note: a size mismatch or empty input throws std::invalid_argument. To evaluate many points at the same sample set, barycentricLagrangeWeights / barycentricLagrangeEvaluate, which reuse the weights, are faster. Beware Runge's phenomenon at high degree with equally spaced points (see floaterHormannWeights).

// Example: evaluate the parabola through (0,0),(1,1),(2,4) at x=1.5
std::vector<double> xs{0.0, 1.0, 2.0}, ys{0.0, 1.0, 4.0};
double v = lagrange_interpolation<double>(xs, ys, 1.5);
// Result: v = 2.25  (= 1.5^2)

newton_interpolation_coefficients

template<concepts::Field T>
std::vector<T> newton_interpolation_coefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

Behaviour: builds the Newton forward divided-difference table and returns the coefficients $\{c_i\}$ of the interpolating polynomial $p(x) = c_0 + c_1 (x - x_0) + c_2 (x - x_0)(x - x_1) + \cdots$. The coefficients are passed to newton_interpolation_evaluate for evaluation.

Parameters:

Use when: the sample points are fixed and the interpolant is evaluated many times. Once the coefficients are computed, each evaluation is $O(n)$ via Horner's method.

newton_interpolation_evaluate

template<concepts::Field T>
T newton_interpolation_evaluate(
    const std::vector<T>& x,
    const std::vector<T>& coeffs,
    T xi);

Behaviour: evaluates the divided-difference coefficients returned by newton_interpolation_coefficients at the point $x_i$ using Horner's method.

Parameters:

std::vector<double> xs{0.0, 1.0, 2.0}, ys{0.0, 1.0, 4.0};
auto c = newton_interpolation_coefficients<double>(xs, ys);
double v = newton_interpolation_evaluate<double>(xs, c, 1.5);
// Result: v = 2.25  (the same polynomial as Lagrange interpolation)

linear_interpolate

template<concepts::Field T>
T linear_interpolate(T x0, T y0, T x1, T y1, T xi);

Behaviour: returns the value at $x_i$ on the line through the two points $(x_0, y_0)$, $(x_1, y_1)$ ($y_0 + t(y_1 - y_0)$ with $t = (x_i - x_0)/(x_1 - x_0)$). When $x_0 \approx x_1$, it returns the average $(y_0 + y_1)/2$ of the two points to avoid the degenerate case.

Parameters:

bilinear_interpolate

template<concepts::Field T>
T bilinear_interpolate(
    T x, T y,
    T x1, T y1, T x2, T y2,
    T q11, T q12, T q21, T q22);

Behaviour: bilinearly interpolates the value at the point $(x, y)$ from the four corner values of the rectangle $[x_1, x_2] \times [y_1, y_2]$. It interpolates linearly along the bottom and top edges, then linearly again in the vertical direction.

Parameters:

Note: if the evaluation point lies outside the rectangle $[x_1, x_2] \times [y_1, y_2]$, it throws std::invalid_argument.

stepInterpolation

template<concepts::OrderedField T>
T stepInterpolation(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T xi);

Behaviour: step interpolation (zero-order hold). Returns the left-endpoint value $y_i$ of the interval $[x_i, x_{i+1})$ containing $x_i$ unchanged. The interval lookup is a binary search, $O(\log n)$.

Parameters:

Note: out of range, it returns the value of the nearest endpoint ($y_0$ for $x_i \le x_0$, $y_{n-1}$ for $x_i \ge x_{n-1}$).

barycentricLagrangeWeights

template<concepts::Field T>
std::vector<T> barycentricLagrangeWeights(
    const std::vector<T>& x);

Behaviour: computes the barycentric Lagrange interpolation weights $w_j = 1 / \prod_{i \ne j} (x_j - x_i)$. The weights depend only on the sample points, not on the values $y$, so they can be computed once and reused.

Parameters:

Use when: plain Lagrange interpolation costs $O(n^2)$ every time, whereas the barycentric form costs $O(n^2)$ for the weights plus $O(n)$ per evaluation. Advantageous when evaluating many points at the same sample set.

barycentricLagrangeEvaluate

template<concepts::Field T>
T barycentricLagrangeEvaluate(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& weights,
    T xi);

Behaviour: evaluates the interpolated value with the second barycentric formula $p(x_i) = \dfrac{\sum_j w_j y_j / (x_i - x_j)}{\sum_j w_j / (x_i - x_j)}$. If $x_i$ coincides with one of the sample points, that $y_j$ is returned directly.

Parameters:

std::vector<double> xs{0.0, 1.0, 2.0}, ys{0.0, 1.0, 4.0};
auto w = barycentricLagrangeWeights<double>(xs);
double v = barycentricLagrangeEvaluate<double>(xs, ys, w, 1.5);
// Result: v = 2.25

Cubic Spline

A cubic spline is a piecewise cubic curve that joins adjacent intervals smoothly with $C^2$ continuity. sangi offers two entry points.

• cubic_spline_coefficients + cubic_spline_evaluate, which return a coefficient array (Interpolation.hpp)
• cubicSpline, which returns a result struct bundling the coefficients and an eval method (BSpline.hpp)

SplineBoundaryCondition

enum class SplineBoundaryCondition {
    Natural,   // natural spline (second derivative is 0 at both ends)
    Clamped    // clamped spline (first-derivative value specified at both ends)
};

cubic_spline_coefficients

template<concepts::Field T>
std::vector<std::array<T, 4>> cubic_spline_coefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    int boundary_condition = 0,
    const std::pair<T, T>& bc_values = { 0, 0 });

Behaviour: computes the coefficients of the piecewise cubic spline through the sample points. The coefficients of each interval are $\{a, b, c, d\}$, so on interval $i$ we have $S_i(x) = a + b(x - x_i) + c(x - x_i)^2 + d(x - x_i)^3$. Internally it solves a tridiagonal system with the Thomas algorithm (forward elimination + back substitution). It returns $n-1$ std::array<T, 4> (one per interval).

Parameters:

Note: the span overload takes the third argument as a SplineBoundaryCondition enumerator. A size mismatch, fewer than 2 points, or non-increasing $x$ throws std::invalid_argument. The returned coefficients are evaluated with cubic_spline_evaluate.

cubic_spline_evaluate

template<concepts::Field T>
T cubic_spline_evaluate(
    const std::vector<T>& x,
    const std::vector<std::array<T, 4>>& coeffs,
    T xi);

Behaviour: evaluates the spline coefficients at the point $x_i$. It locates the interval containing $x_i$ by linear search and returns $S_i(x_i) = a + b\,\delta + c\,\delta^2 + d\,\delta^3$ ($\delta = x_i - x_i^{\text{left}}$). Out of range, it returns the constant term of the last interval.

Parameters:

Use when: not only cubic_spline_coefficients, but also pchipCoefficients, catmullRomCoefficients, modifiedAkimaCoefficients, and smoothInterpolationCoefficients all return coefficients in the same $\{a, b, c, d\}$ form, so this function can evaluate any of them.

std::vector<double> xs{0.0, 1.0, 2.0, 3.0}, ys{0.0, 1.0, 0.0, 1.0};
auto coeffs = cubic_spline_coefficients<double>(xs, ys);  // default: natural boundary
double v = cubic_spline_evaluate<double>(xs, coeffs, 1.5);
// Result: v ~ around 0.5 (passing smoothly between sample points (1,1)-(2,0))

cubicSpline

template<typename T>
struct CubicSplineResult {
    std::vector<T> a, b, c, d;  // per-interval coefficients: S_i(x) = a + b(x-x_i) + c(x-x_i)^2 + d(x-x_i)^3
    std::vector<T> x;           // knots
    T eval(T t) const;          // evaluate the interpolated value
};

template<typename T>
CubicSplineResult<T> cubicSpline(
    const std::vector<T>& xs,
    const std::vector<T>& ys);

Behaviour: the natural cubic spline interpolation ($S'' = 0$ at both ends) provided by BSpline.hpp. It returns a CubicSplineResult<T> bundling the coefficients and the evaluation method eval, so you do not have to carry the coefficient array and $x$ around separately.

Parameters:

Returns: CubicSplineResult<T>. Use result.eval(t) to get the interpolated value at any point.

std::vector<double> xs{0.0, 1.0, 2.0, 3.0}, ys{0.0, 1.0, 4.0, 9.0};
auto spline = cubicSpline<double>(xs, ys);
double v = spline.eval(1.5);
// Result: v ~ around 2.3 (passing smoothly between knots (1,1)-(2,4))

Shape-Preserving & Local Interpolation

The following functions return per-interval cubic (and, in one case, quintic) coefficients, evaluated with cubic_spline_evaluate (quintic with quinticHermiteEvaluate). Because they fix slopes locally rather than performing a global optimization, they tend to suppress oscillation (overshoot) better than the standard cubic spline.

Function reference

pchipCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> pchipCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

Behaviour: PCHIP (Piecewise Cubic Hermite Interpolating Polynomial). The Fritsch-Carlson algorithm fixes the slope at each knot to preserve monotonicity on monotone intervals and suppress overshoot. Interior slopes are estimated with the harmonic mean and endpoint slopes with Bessel's one-sided three-point formula, then limited to satisfy the monotonicity condition $\alpha^2 + \beta^2 \le 9$. It returns the per-interval $\{a, b, c, d\}$, evaluated with cubic_spline_evaluate.

Parameters:

Use when: you want to interpolate monotone data (cumulative distributions, concentration profiles, etc.) smoothly but without non-physical under-/overshoot. Fewer than 2 points or a size mismatch returns an empty vector.

std::vector<double> xs{0.0, 1.0, 2.0, 3.0}, ys{0.0, 0.0, 1.0, 1.0};  // step-like monotone data
auto coeffs = pchipCoefficients<double>(xs, ys);
double v = cubic_spline_evaluate<double>(xs, coeffs, 1.5);
// Result: 0 <= v <= 1 (monotonicity preservation means it never overshoots the range)

modifiedAkimaCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> modifiedAkimaCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

Behaviour: a version of Akima's local slope estimate with the Boost.Math-compatible modification. By adding a tiny quantity $\varepsilon$ to the absolute value of the slope difference $|\delta_{i+1} - \delta_i|$, it avoids division by zero and improves stability on flat regions and equally spaced data. It returns the per-interval $\{a, b, c, d\}$, evaluated with cubic_spline_evaluate.

Parameters:

Use when: you want to confine the influence of outliers locally. A single anomalous value is unlikely to propagate to distant intervals (the hallmark of Akima). Fewer than 2 points or a size mismatch returns an empty vector.

catmullRomCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> catmullRomCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

Behaviour: the Catmull-Rom spline. A $C^1$-continuous interpolating spline whose tangent at each knot is found by a difference of the two neighboring knots. The endpoints use one-sided differences. It returns the per-interval $\{a, b, c, d\}$, evaluated with cubic_spline_evaluate.

Parameters:

Use when: keyframe interpolation and curve drawing in computer graphics. Fewer than 2 points or a size mismatch returns an empty vector.

smoothInterpolationCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 4>> smoothInterpolationCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y);

Behaviour: a piecewise interpolation that blends the slopes of adjacent intervals to smooth the joins. The blending formula $A_{\text{mix}} = (|A_0| A_1 + |A_1| A_0)/(|A_0| + |A_1|)$ reproduces a straight line when three points are collinear, yields a blend of 0 when one side is flat, and yields a blend of 0 at a peak (slopes of opposite sign). The first and last intervals are quadratic with the endpoint slope constrained; the interior intervals are cubic with both end slopes constrained. It returns the per-interval $\{a, b, c, d\}$, evaluated with cubic_spline_evaluate.

Parameters:

Note: if $x$ is not strictly increasing, it throws std::invalid_argument.

floaterHormannWeights

template<concepts::OrderedField T>
std::vector<T> floaterHormannWeights(
    const std::vector<T>& x,
    std::size_t d = 3);

Behaviour: computes the weights for Floater-Hormann barycentric rational interpolation. This is a rational interpolant blending local polynomials of degree $d$ that is stable and free of Runge's phenomenon even on equally spaced points. $d = 0$ corresponds to piecewise constant, and $d = n-1$ to polynomial interpolation.

Parameters:

Use when: you want to interpolate high-degree equally spaced data with a polynomial-like fit but avoid Runge's phenomenon. The weights are passed to floaterHormannEvaluate.

floaterHormannEvaluate

template<concepts::OrderedField T>
T floaterHormannEvaluate(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& weights,
    T xi);

Behaviour: evaluates the Floater-Hormann rational interpolation with the barycentric formula $r(x_i) = \dfrac{\sum_k w_k y_k / (x_i - x_k)}{\sum_k w_k / (x_i - x_k)}$. If the evaluation point coincides with (or is extremely close to) a knot, that $y_k$ is returned.

Parameters:

quinticHermiteCoefficients

template<concepts::OrderedField T>
std::vector<std::array<T, 6>> quinticHermiteCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& dy,
    const std::vector<T>& d2y);

Behaviour: quintic Hermite interpolation specifying the value $f$, first derivative $f'$, and second derivative $f''$ at each knot. $C^2$-continuous. It returns 6 coefficients $\{a, b, c, d, e, f\}$ per interval, evaluated as $S(t) = a + b\,t + c\,t^2 + d\,t^3 + e\,t^4 + f\,t^5$ ($t = x - x_i$).

Parameters:

Use when: interpolating trajectories or motion profiles where you want to specify the curvature as well. A size mismatch or fewer than 2 points returns an empty vector. Evaluate with quinticHermiteEvaluate (since the coefficient length is 6, cubic_spline_evaluate cannot be used).

quinticHermiteEvaluate

template<concepts::OrderedField T>
T quinticHermiteEvaluate(
    const std::vector<T>& x,
    const std::vector<std::array<T, 6>>& coeffs,
    T xi);

Behaviour: evaluates the coefficients from quinticHermiteCoefficients at the point $x_i$. It locates the interval by binary search and clamps to the end interval when out of range. The quintic polynomial is evaluated by Horner's method.

Parameters:

B-spline

A B-spline represents a curve as a linear combination of basis functions with compact support. sangi provides two families of basis evaluation.

• bsplineBasis(degree, i, x, knots) from Interpolation.hpp (Cox-de Boor recursion, OrderedField), bsplineCoefficients, bsplineEvaluate
• bsplineBasis(i, p, t, knots) from BSpline.hpp (De Boor recursion, int degree), bsplineBasisAll, uniformKnots, deBoor, bsplineRegression

Both headers contain a function named bsplineBasis, but the argument orders differ (the former takes a std::span<const T> knot vector and a std::size_t degree; the latter takes a std::vector<T> knot vector and an int degree).

Function reference

bsplineBasis (Interpolation.hpp)

template<concepts::OrderedField T>
T bsplineBasis(
    std::size_t degree,
    std::size_t i,
    T x,
    std::span<const T> knots);

Behaviour: evaluates the B-spline basis function $B_{i,\text{degree}}(x)$ by Cox-de Boor recursion. At degree 0 it is the indicator function of the interval $[\text{knot}_i, \text{knot}_{i+1})$ (only the last interval includes its right endpoint).

Parameters:

bsplineCoefficients

template<concepts::OrderedField T>
std::pair<std::vector<T>, std::vector<T>> bsplineCoefficients(
    const std::vector<T>& x,
    const std::vector<T>& y,
    std::size_t degree = 3);

Behaviour: returns the pair {knots, coefficients} of the knot vector and control-point coefficients of the B-spline that interpolates the sample points. It places clamped knots of multiplicity $\text{degree}+1$ at both ends and chooses the interior knots as averaging knots. It assembles the collocation matrix $B_{j,\text{degree}}(x_i)$ and solves the linear system by LU decomposition.

Parameters:

Note: a size mismatch, or fewer than $\text{degree}+1$ points, throws std::invalid_argument. The returned knots and coefficients are passed to bsplineEvaluate.

bsplineEvaluate

template<concepts::OrderedField T>
T bsplineEvaluate(
    const std::vector<T>& knots,
    const std::vector<T>& coefficients,
    std::size_t degree,
    T x);

Behaviour: evaluates the B-spline interpolated value as the linear combination of control-point coefficients and bases $\sum_i c_i B_{i,\text{degree}}(x)$.

Parameters:

std::vector<double> xs{0.0, 1.0, 2.0, 3.0, 4.0}, ys{0.0, 1.0, 0.0, 1.0, 0.0};
auto [knots, coeffs] = bsplineCoefficients<double>(xs, ys, 3);
double v = bsplineEvaluate<double>(knots, coeffs, 3, 2.5);
// value on the cubic B-spline through the sample points (between xs[2]=2 and xs[3]=3)

bsplineBasis (BSpline.hpp)

template<typename T>
T bsplineBasis(size_t i, int p, T t, const std::vector<T>& knots);

Behaviour: evaluates the basis $N_{i,p}(t)$ by De Boor recursion (the same Cox-de Boor recurrence as the Interpolation.hpp version, but a separate overload taking the degree as int and the knots as std::vector).

Parameters:

bsplineBasisAll

template<typename T>
std::vector<T> bsplineBasisAll(int p, T t, const std::vector<T>& knots);

Behaviour: evaluates all bases $N_{0,p}(t), \ldots, N_{n-1,p}(t)$ ($n = \text{knots.size}() - p - 1$) at the point $t$ in one pass. When $t$ is at the right end of the knot vector, only the last basis is set to 1 (to avoid missing the right endpoint).

Parameters:

uniformKnots

template<typename T>
std::vector<T> uniformKnots(int n, int p, T a = T{0}, T b = T{1});

Behaviour: generates a clamped uniform knot vector over the interval $[a, b]$ (multiplicity $p+1$ at both ends, equally spaced in the interior). The total number of knots is $n + p + 1$.

Parameters:

deBoor

template<typename T>
std::vector<T> deBoor(
    const std::vector<std::vector<T>>& controlPoints,
    const std::vector<T>& knots,
    int p,
    T t);

Behaviour: evaluates a single point on the B-spline curve $C(t) = \sum_i N_{i,p}(t) P_i$ by de Boor's algorithm. It locates the knot span and recursively contracts the triangular table to return a point on the curve (a vector of dimension $d$).

Parameters:

Returns: a point on the curve (a std::vector<T> of the same dimension $d$ as the control points).

bsplineRegression

template<typename T>
struct BSplineRegressionResult {
    std::vector<T> coefficients;  // B-spline coefficients
    std::vector<T> knots;         // knot vector
    int degree;                   // degree
    T eval(T t) const;            // evaluate the predicted value
};

template<typename T>
BSplineRegressionResult<T> bsplineRegression(
    const std::vector<T>& x,
    const std::vector<T>& y,
    int nBasis = 10,
    int degree = 3);

Behaviour: least-squares regression onto $n_{\text{basis}}$ B-spline bases. It builds the design matrix $B$ with uniform knots and solves the normal equations $B^\top B\, c = B^\top y$ (with a tiny regularization added) by Gaussian elimination. This is a smooth fit, not an interpolation (it does not pass through every point). Use the result type's eval(t) to get the predicted value.

Parameters:

Use when: you want to draw a smooth trend curve through noisy data. Adjust the balance between smoothness and fit with the number of bases nBasis.

Trigonometric Interpolation

cardinalTrigonometricInterpolate

template<concepts::OrderedField T>
T cardinalTrigonometricInterpolate(
    const std::vector<T>& y,
    T period,
    T xi);

Behaviour: interpolates equally spaced periodic data $y[0..N-1]$ (sample points $x_k = k \cdot \text{period}/N$) with a trigonometric polynomial. It computes the Dirichlet-kernel approach $S(x) = \frac{1}{N}\sum_k y_k\, D_N(x - x_k)$ directly (no FFT, suited to small to medium data). It switches kernels depending on whether $N$ is even or odd.

Parameters:

Use when: reconstructing periodic data that underlies periodic signals and spectral methods. It reproduces the original values exactly at the sample points.

// Example: sample a sine wave of period 2pi at 8 points and interpolate
std::vector<double> y(8);
for (int k = 0; k < 8; ++k) y[k] = std::sin(2.0 * M_PI * k / 8.0);
double v = cardinalTrigonometricInterpolate<double>(y, 2.0 * M_PI, 0.3);
// Result: v ~ sin(0.3) (reconstructed smoothly with a trigonometric polynomial between samples)

Spline Fitting (Smoothing)

Rather than an interpolation that passes exactly through the sample points, this is a least-squares fit that trades off "fit" against "smoothness" via a smoothing parameter $s$. It corresponds to FITPACK's (Dierckx's) splrep/splev. There are three kinds: 1D, parametric curve, and 2D surface.

Function reference

SplineFitResult / splineFit

template<concepts::OrderedField T>
struct SplineFitResult {
    std::vector<T> knots;        // knot vector
    std::vector<T> coefficients; // B-spline coefficients
    std::size_t degree;          // spline degree
    T smoothing;                 // actual smoothing parameter
    T residual;                  // residual sum of squares
};

template<concepts::OrderedField T>
SplineFitResult<T> splineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& w = {},
    T s = T(-1),
    std::size_t degree = 3);

Behaviour: a smoothing B-spline fit that minimizes $\sum_i w_i (y_i - S(x_i))^2 + s \int (S''(x))^2 dx$ for the data $(x_i, y_i)$. $s = 0$ gives interpolation (passing through every point), $s > 0$ gives smoothing (noise suppression). It automatically selects the number of knots according to the smoothing parameter.

Parameters:

Note: fewer than $\text{degree}+1$ points, a size mismatch, or an invalid weight length throws std::invalid_argument. The result is evaluated with splineEval.

splineEval

template<concepts::OrderedField T>
T splineEval(const SplineFitResult<T>& fit, T x);

template<concepts::OrderedField T>
std::vector<T> splineEval(const SplineFitResult<T>& fit, std::span<const T> xs);

Behaviour: evaluates the result of splineFit at arbitrary points (equivalent to FITPACK's splev). There is a single-point overload and a multi-point overload that evaluates many points at once.

Parameters:

// Smoothing fit of noisy data
std::vector<double> xs = /* ascending x */;
std::vector<double> ys = /* noisy y */;
auto fit = splineFit<double>(xs, ys);          // automatic s, cubic
double v = splineEval<double>(fit, 1.5);        // value on the smoothing curve
// fit.residual holds the residual sum of squares, fit.smoothing the s that was used

ParametricSplineResult / parametricSplineFit

template<concepts::OrderedField T>
struct ParametricSplineResult {
    SplineFitResult<T> x_fit;   // fit of x(t)
    SplineFitResult<T> y_fit;   // fit of y(t)
    std::vector<T> t;            // parameter values
};

template<concepts::OrderedField T>
ParametricSplineResult<T> parametricSplineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    T s = T(-1),
    std::size_t degree = 3);

Behaviour: fits a planar curve $(x(t), y(t))$ using an arc-length parameter (cumulative chord length normalized to $[0, 1]$). It fits $x(t)$ and $y(t)$ separately, each with splineFit. $x$ need not be monotone (closed curves / loops are allowed).

Parameters:

Note: fewer than 2 points or a size mismatch throws std::invalid_argument. Evaluate with parametricSplineEval.

parametricSplineEval

template<concepts::OrderedField T>
std::pair<T, T> parametricSplineEval(
    const ParametricSplineResult<T>& fit,
    T t);

Behaviour: returns the point $(x(t), y(t))$ on the curve at parameter $t \in [0, 1]$ as a std::pair<T, T>.

Parameters:

SurfaceSplineResult / surfaceSplineFit

template<concepts::OrderedField T>
struct SurfaceSplineResult {
    std::vector<T> knots_x;      // knot vector in the x direction
    std::vector<T> knots_y;      // knot vector in the y direction
    std::vector<T> coefficients; // coefficients (nx x ny row-major flat array)
    std::size_t nx;             // number of bases in the x direction
    std::size_t ny;             // number of bases in the y direction
    std::size_t degree;         // spline degree
};

template<concepts::OrderedField T>
SurfaceSplineResult<T> surfaceSplineFit(
    const std::vector<T>& x,
    const std::vector<T>& y,
    const std::vector<T>& z,
    std::size_t degree = 3);

Behaviour: fits grid data $z[i][j] = f(x_i, y_j)$ with a tensor-product B-spline $S(x, y) = \sum_i \sum_j c_{ij} B_i(x) B_j(y)$. A two-stage decomposition: it first fits each $x$ row in the $y$ direction, then fits those coefficients in the $x$ direction.

Parameters:

Note: if the number of points in either direction is fewer than $\text{degree}+1$, or if the size of $z$ is not $m_x m_y$, it throws std::invalid_argument. Evaluate with surfaceSplineEval.

surfaceSplineEval

template<concepts::OrderedField T>
T surfaceSplineEval(const SurfaceSplineResult<T>& fit, T x, T y);

Behaviour: evaluates the tensor-product spline $S(x, y) = \sum_i \sum_j c_{ij} B_i(x) B_j(y)$ at the point $(x, y)$.

Parameters:

Image Interpolation (2D Resampling)

A family of functions that interpolate a row-major 2D grid (image) at a continuous coordinate $(x, y)$. $x$ is the column-direction and $y$ the row-direction 0-based coordinate. Boundaries are clamped.

Function reference

nearestNeighbor2D

template<concepts::OrderedField T>
T nearestNeighbor2D(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y);

Behaviour: returns the value of the nearest grid point unchanged (rounds the coordinates to an index, clamping the range). Used for downscaling images or upscaling pixel art.

Parameters:

bicubicInterpolate

template<concepts::OrderedField T>
T bicubicInterpolate(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y);

Behaviour: smoothly interpolates by weighting the $4 \times 4 = 16$ neighboring grid points with the Catmull-Rom kernel ($a = -0.5$). It corresponds to OpenCV's INTER_CUBIC and Photoshop's bicubic method. Boundaries clamp the index.

Parameters:

Use when: you want higher quality than nearest-neighbor or bilinear when enlarging or shrinking photos.

lanczosInterpolate2D

template<concepts::OrderedField T>
T lanczosInterpolate2D(
    std::span<const T> grid,
    std::size_t rows, std::size_t cols,
    T x, T y,
    int a = 3);

Behaviour: high-quality resampling based on sinc. It applies the Lanczos-$a$ kernel $L(x) = \mathrm{sinc}(x)\,\mathrm{sinc}(x/a)$ ($|x| < a$, 0 otherwise) over a $2a \times 2a$ window. It normalizes so that the weight sum does not break down at boundaries. The weights are computed in the $T$ type to avoid loss of precision (preserving precision even with the arbitrary-precision Float).

Parameters:

Use when: high-quality image down-/upscaling. It balances ringing and sharpness better than bicubic and is a staple for photo resampling.

Example

#include <math/interpolation/Interpolation.hpp>
#include <math/interpolation/BSpline.hpp>
#include <iostream>
#include <vector>
using namespace sangi;

int main() {
    // ---- 1) interpolate a point sequence with cubicSpline and eval ----
    std::vector<double> xs{0.0, 1.0, 2.0, 3.0, 4.0};
    std::vector<double> ys{0.0, 1.0, 4.0, 9.0, 16.0};   // samples of y = x^2
    auto spline = cubicSpline<double>(xs, ys);
    std::cout << "spline(2.5) = " << spline.eval(2.5) << '\n';
    // Output: spline(2.5) ~ around 6.25 (2.5^2 = 6.25)

    // ---- 2) cubic_spline_coefficients / evaluate (coefficient-array version) ----
    auto coeffs = cubic_spline_coefficients<double>(xs, ys);  // natural boundary
    std::cout << "eval(1.5) = "
              << cubic_spline_evaluate<double>(xs, coeffs, 1.5) << '\n';

    // ---- 3) smoothing with splineFit (noise suppression) ----
    std::vector<double> xn{0,1,2,3,4,5,6,7,8,9};
    std::vector<double> yn{0.1,0.9,2.1,2.9,4.2,4.8,6.1,6.9,8.2,8.9};  // y ~ x + noise
    auto fit = splineFit<double>(xn, yn);   // automatic s, cubic
    std::cout << "fit(4.5) = " << splineEval<double>(fit, 4.5) << '\n';
    std::cout << "residual = " << fit.residual << '\n';

    // ---- 4) monotonicity-preserving interpolation with PCHIP ----
    std::vector<double> xm{0,1,2,3}, ym{0,0,1,1};   // step-like monotone data
    auto pc = pchipCoefficients<double>(xm, ym);
    std::cout << "pchip(1.5) = "
              << cubic_spline_evaluate<double>(xm, pc, 1.5) << '\n';
    // monotonicity preservation keeps it within 0 to 1 (no overshoot)
}

Related Mathematical Background

The following articles explain the mathematical concepts underlying the Interpolation module.

• Interpolation — Reconstructing intermediate values from sample points
• Lagrange Interpolating Polynomial — Polynomial interpolation and the Runge phenomenon
• Cubic Splines — $C^2$-continuous piecewise cubics and the tridiagonal system
• B-splines — Cox-de Boor basis and locality
• Hermite Interpolation — Matching derivatives as well as values
• Least Squares — Foundation of spline fitting and regression