Interpolation: B-splines and Smoothing Spline Fitting

Overview

A B-spline (basis spline) is a framework that expresses a piecewise polynomial—low-degree polynomials joined smoothly at knots—as a linear combination of basis functions $N_{i,p}(t)$ with local support. The curve can be written as a weighted sum of control points $P_i$: $C(t) = \sum_i N_{i,p}(t)\, P_i$.

There are three main advantages over polynomial interpolation:

• Locality: each basis $N_{i,p}$ is non-zero on only $p+1$ knot spans. Moving a single control point changes only the nearby part of the curve and does not propagate globally.
• Convex hull property: the basis functions are non-negative and sum to $1$ (partition of unity), so at every point the curve lies within the convex hull of the control points. This keeps the shape from misbehaving and makes it numerically stable as well.
• Arbitrary degree: you can choose the degree $p$. Low-degree piecewise polynomials avoid the Runge phenomenon (violent oscillation near the ends) that high-degree polynomial interpolation produces on equally spaced data.

Related API: bsplineCoefficients, bsplineEvaluate, bsplineRegression, splineFit.

Cox-de Boor Basis

The B-spline basis functions of degree $p$ are defined recursively over a knot vector $\{t_0 \le t_1 \le \dots \le t_m\}$ by the Cox-de Boor recurrence. It starts from degree 0 (piecewise constant):

$$N_{i,0}(t) = \begin{cases} 1 & t_i \le t < t_{i+1} \\ 0 & \text{otherwise} \end{cases}$$

For degree $p \ge 1$, two adjacent lower-degree bases are blended with linear weights in $t$:

$$N_{i,p}(t) = \frac{t - t_i}{t_{i+p} - t_i}\, N_{i,p-1}(t) + \frac{t_{i+p+1} - t}{t_{i+p+1} - t_{i+1}}\, N_{i+1,p-1}(t)$$

A term whose denominator becomes $0$ (a repeated knot) is handled by the convention that the entire term is taken as $0$. sangi's bsplineBasis implements this recurrence directly and drops a term when its denominator is at or below machine epsilon.

Support interval

$N_{i,p}(t)$ is identically $0$ outside the interval $[t_i,\ t_{i+p+1})$. That is, each basis affects at most $p+1$ knot spans. This is the basis of locality, and at evaluation time only the non-zero bases need to be considered, which keeps the cost down.

Clamped uniform knot vector

To make the curve pass reliably through the endpoints, a clamped knot vector is used, with the knots at both ends repeated with multiplicity $p+1$. On an interval $[a, b]$ with $n$ basis functions, the total number of knots is $n + p + 1$:

$$\underbrace{a, \dots, a}_{p+1},\ \underbrace{t_{p+1}, \dots, t_{n-1}}_{\text{interior knots}},\ \underbrace{b, \dots, b}_{p+1}$$

sangi's uniformKnots places the interior knots at equal spacing. For interpolation (bsplineCoefficients), it also offers a scheme that places the interior knots at the averaged positions of the data points (averaged knots) to improve the conditioning of the collocation matrix.

De Boor Evaluation

Given control points $P_i$ and a knot vector, the curve value at parameter $t$ is defined as

$$C(t) = \sum_{i} N_{i,p}(t)\, P_i$$

The De Boor algorithm computes this without explicitly evaluating the basis functions. First it locates the knot span $[t_k,\ t_{k+1})$ containing $t$, and takes the relevant $p+1$ control points $P_{k-p}, \dots, P_k$ as initial values, repeating interpolation in a triangular scheme.

At each stage $r = 1, \dots, p$, using the coefficient

$$\alpha = \frac{t - t_{k+1+j-p-1}}{t_{k+1+j-r} - t_{k+1+j-p-1}}$$

it updates

$$d_j^{(r)} = (1 - \alpha)\, d_{j-1}^{(r-1)} + \alpha\, d_j^{(r-1)}$$

and the final stage value $d_p^{(p)}$ is the curve value $C(t)$.

Numerical stability

Each update is a convex combination with coefficient $\alpha \in [0, 1]$. Because the value is always interpolated as an internal division of two points, coefficients are not amplified and rounding error is unlikely to accumulate. De Boor's triangular algorithm is numerically more robust than the naive method of recursively expanding the bases and summing. sangi's deBoor sets $\alpha = 0$ when the denominator is at or below machine epsilon, so it does not break down even with repeated knots.

B-spline Interpolation and Regression

Interpolation

Find coefficients $c_j$ that pass exactly through all data points $(x_i, y_i)$ ($i = 0, \dots, n-1$). With the basis matrix (collocation matrix) $B_{ij} = N_{j,p}(x_i)$, the interpolation condition $\sum_j c_j N_{j,p}(x_i) = y_i$ becomes the square linear system

$$B\, c = y$$

By the locality of the basis, $B$ is banded and can be solved efficiently. sangi's bsplineCoefficients solves this system with clamped knots via LU decomposition. Furthermore, the endpoint rows ($x_0$ and $x_{n-1}$) are set explicitly using the facts $N_{0,p}(x_0) = 1$ and $N_{n-1,p}(x_{n-1}) = 1$, preventing numerical error at the endpoints.

Regression

If the number of basis functions $M$ is taken smaller than the number of data points $n$, $B$ becomes an $n \times M$ tall matrix and an exact solution generally does not exist. So the least-squares solution that minimizes the sum of squared residuals $\|y - Bc\|^2$ is found from the normal equations

$$B^\top B\, c = B^\top y$$

where $B^\top B$ is an $M \times M$ symmetric positive semi-definite matrix. The fewer the basis functions, the smoother the curve and the more the data's noise is averaged out.

sangi's bsplineRegression constructs the normal equations and then adds a small regularization term on the diagonal,

$$(B^\top B + \lambda I)\, c = B^\top y, \qquad \lambda = 10^{-10}$$

solving it by Gaussian elimination. $\lambda$ is a ridge term that prevents numerical breakdown when $B^\top B$ is nearly singular (for example, when some basis is not excited over part of the data range), and it is taken small enough to have almost no effect on the shape of the fit.

Smoothing fit and the smoothing parameter $s$

A smoothing spline continuously bridges between interpolation ($s=0$) and strong smoothing. It constrains the weighted sum of squared residuals, with weights $w_i$, by the smoothing parameter $s$ as an upper bound:

$$\sum_{i} w_i\, \bigl(y_i - S(x_i)\bigr)^2 \le s$$

• $s = 0$: interpolation passing exactly through all data points.
• large $s$: tolerate the residual, reduce the knots, and get a smooth curve. The effect of noise is suppressed.

$s$ governs the trade-off between residual (fidelity to the data) and smoothness. sangi's splineFit increases the interior knots a little at a time, performs a weighted least-squares fit at each pass, and stops once the residual drops to or below $s$ (to avoid overfitting, it terminates early if the residual starts to worsen). If $s$ is specified as negative, it follows Dierckx's recommendation and auto-selects $s = n$ (the number of data points) as the default.

The least squares at each candidate knot count is solved as the weighted normal equations using the weight matrix $W = \mathrm{diag}(w_i)$:

$$(B^\top W B)\, c = B^\top W y$$

The implementation forms $W^{1/2}B$ and $W^{1/2}y$, then assembles the normal equations and solves them by LU decomposition.

Parametric Curve and Surface Splines

Chord-length parametric curve

A closed curve or a folded trajectory where $y$ is not a function of $x$ cannot be represented in the form $y = f(x)$. In this case, $x(t)$ and $y(t)$ are fit by separate splines against a common parameter $t$.

Chord-length parameterization is used for the parameter. The Euclidean distances between adjacent data points are accumulated and normalized by the total length:

$$t_0 = 0, \qquad t_i = t_{i-1} + \sqrt{(x_i - x_{i-1})^2 + (y_i - y_{i-1})^2}, \qquad t_i \leftarrow t_i / t_{n-1}$$

This way the parameter advances slowly over densely sampled intervals, giving a fit with little bias regardless of the curve's shape. sangi's parametricSplineFit fits $x(t)$ and $y(t)$ separately with splineFit, and evaluation returns the pair of their values $(x(t), y(t))$.

2D tensor-product surface spline

To represent grid data $z_{ij} = f(x_i, y_j)$ with a smooth surface, a product of B-spline bases in the two directions (a tensor product) is used:

$$S(x, y) = \sum_i \sum_j c_{ij}\, N_i(x)\, N_j(y)$$

The coefficients $c_{ij}$ are obtained by a two-stage 1D solve. First, for each $x$ row, the B-spline coefficients in the $y$ direction are found (intermediate coefficients $\alpha_{ij}$); then, for each of those $y$ columns, the B-spline coefficients in the $x$ direction are found. A clamped knot vector is built for each of the $x$ and $y$ directions, and each 1D system is solved by LU decomposition. sangi's surfaceSplineFit implements this separated solve, and the evaluator surfaceSplineEval computes the double sum by looping only over the non-zero $x$-direction bases.

Image Interpolation

Image resampling (scaling, shrinking, rotation) is the problem of reconstructing the pixel values $g[r][c]$ on an integer grid at non-integer evaluation points $(x, y)$. The reconstruction is expressed as a convolution with an interpolation kernel $K$, $\hat g(x, y) = \sum_{m, n} K(x - n)\, K(y - m)\, g[m][n]$, and the width and shape of the kernel's support determine the reconstruction quality. sangi provides the following three.

Nearest neighbor

Round the evaluation point to the closest grid point and return that pixel value as is (nearestNeighbor2D). The support is one pixel. It is the fastest, but produces stair-step discontinuities, so it suits uses where blurring is to be avoided, such as integer-factor upscaling of pixel art.

Bicubic (Catmull-Rom kernel)

Weight the $4 \times 4 = 16$ neighboring pixels with a cubic polynomial kernel (bicubicInterpolate). sangi adopts the Catmull-Rom kernel ($a = -1/2$):

$$K(t) = \begin{cases} \tfrac{3}{2}|t|^3 - \tfrac{5}{2}|t|^2 + 1 & |t| \le 1 \\[4pt] -\tfrac{1}{2}|t|^3 + \tfrac{5}{2}|t|^2 - 4|t| + 2 & 1 < |t| < 2 \\[4pt] 0 & |t| \ge 2 \end{cases}$$

This kernel has support $[-2, 2]$ and satisfies the interpolation conditions $K(0)=1$, $K(\pm 1)=K(\pm 2)=0$ at the grid points. Because it is smooth and preserves edges relatively well, it is widely used as a general default for resampling.

Lanczos resampling

The Lanczos kernel, a windowed $\mathrm{sinc}$ function, gives the sharpest reconstruction (lanczosInterpolate2D):

$$K(t) = \begin{cases} \mathrm{sinc}(t)\, \mathrm{sinc}(t/a) & |t| < a \\ 0 & |t| \ge a \end{cases}, \qquad \mathrm{sinc}(t) = \frac{\sin(\pi t)}{\pi t}$$

$a$ is the window width: $a = 2$ uses a $4 \times 4$ window, and $a = 3$ (default) uses a $6 \times 6$ window. The wider the support, the closer to an ideal low-pass filter and the higher the quality, but the computation increases and ringing can appear near contours. sangi's implementation normalizes by the sum of the weights to conserve energy near boundaries.

The relationship between the kernel's support and quality can be summarized as a consistent trade-off: the wider the support, the closer the frequency response is to ideal, while computational cost and ringing both increase.

Comparison Tables

Interpolation vs. regression vs. smoothing

Computational complexity

Rough guidance: use interpolation when the data is error-free, regression or smoothing (tuning $s$) when there is noise, a chord-length parametric fit for closed curves or trajectories, a tensor product for gridded surfaces, and choose among nearest neighbor / bicubic / Lanczos for image resampling based on quality and speed.

References

• de Boor, C. (2001). A Practical Guide to Splines. Revised ed. Springer.
• Dierckx, P. (1993). Curve and Surface Fitting with Splines. Oxford University Press.
• Cox, M. G. (1972). "The numerical evaluation of B-splines". Journal of the Institute of Mathematics and its Applications, 10(2), 134–149.
• Piegl, L., & Tiller, W. (1997). The NURBS Book. 2nd ed. Springer.
• Catmull, E., & Rom, R. (1974). "A class of local interpolating splines". In Computer Aided Geometric Design (pp. 317–326). Academic Press.
• Duchon, C. E. (1979). "Lanczos filtering in one and two dimensions". Journal of Applied Meteorology, 18(8), 1016–1022.