Constructing $C^{k-2}$ and $C^{k-1}$ splines on a regular grid through convolution

Question

I am currently trying to figure out how spline interpolation works on a regular grid.

The simplest spline interpolation basis function I can think of is the box function: $f(x) = 1, \, |x|<1$

$$ \phi_0(x) = \beta_0(x)= \begin{cases} 1 & |x| < \frac{1}{2} \\ \frac{1}{2} & |x| = \frac{1}{2} \\ 0 & |x|>\frac{1}{2} \end{cases} $$

Its shifted copies on a regular grid can be used to achieve piecewise-constant interpolation:

$$f(x) = \sum_k f_k \phi_0\left(\frac{x}{h}-k\right)$$

Going one degree higher the Bernstein and Lagrange basis functions are the same:

$$p(x) = p_0 b_{0,1}(x)a + p_1 b_{1,1}(x) = p_0 x + p_1 (1-x)$$

And the spline interpolating basis function is given by shifting $b_{0,1}$ to $[-1,0]$ resulting in $b_{0,1}'(x) = b_{0,1}(x+1) = x+1 = 1-|x|$ and gluing with $b_{1,1}$ yields:

$$ \phi_1(x) = \beta_1(x)= \begin{cases} 1-|x| & |x| < 1 \\ 0 & |x| \geq 1 \end{cases} $$

Then piecewise-linear functions with $C^0$ continuity are reproduced through:

$$f(x) = \sum_k f_k \phi_1\left(\frac{x}{h}-k\right)$$

Starting with degree 2 the Lagrange and Bernstein functions start to differ. But gluing together the Bernstein functions yields:

$$ \beta_2(x) = \begin{cases} \frac{3}{4}-x^2 & |x| < \frac{1}{2} \\ \frac{1}{2}(\frac{3}{2}-|x|)^2 & \frac{1}{2} \leq |x| < \frac{3}{2} \\ 0 & |x| > \frac{3}{2} \end{cases} \quad \beta_3(x) = \begin{cases} \frac{2}{3}-x^2+\frac{1}{2}|x|^3 & |x| < 1 \\ \frac{1}{6}(2-|x|)^3 & 1 \leq |x| < 2 \\ 0 & |x| > 2 \end{cases} $$

And generally it seems like those can be constructed using convolution: $\beta_n = \beta_{n-1} * \beta_0$. My first question is how this gluing together is motivated, and how can this convolution formula be derived from it. The spline basis functions arising from the Bernstein polynomials are not interpolating, but can be made to be by using:

$$f(x) = \sum_k c_k \beta_n\left(\frac{x}{h}-k\right), \quad f_i = \sum_k c_k \beta_n\left(i-k\right)$$

I believe the above is supposed to yield $C^{n-1}$ continuity, what would the proof for this be? Finally, can the same approach be used for gluing the Lagrange functions, in order to construct interpolating basis functions $\phi_n$?

I have found this question: What is the relationship between cubic B-splines and cubic splines?, however it lacked details on why such gluing made sense, also didn't cover the convolution, nor was the Lagrange variant covered.

Answer 1

I have derived the interpolating syntehsis function for $\phi_k$ based on Lagrange interpolation. Unfortunately this interpolation is not $C^{k-2}$, I believe it is only $C^0$, although I haven't verified this formally. I have looked at the resulting synthesis functions for a regular grid and those were $C^0$ in multiple examples.

Interpolant for odd degrees.

For interpolation of odd degree $k$, if we are given a sequence of points $x_1, \ldots, x_n$ then it is similar to the BSpline case. For an interval $[x_i,x_{i+1}]$ one considers the $k+1$ points:

$$x_{i-(k+1)/2+1}, \ldots, x_i, \ldots, x_{i+(k+1)/2},$$

or re-indexing with $a_i = i-(k+1)/2+1$:

$$x_{a_i+0}, \ldots, x_{a_i+(k+1)/2}, \ldots, x_{a_i+k}.$$

Then the interpolant in $[x_i, x_{i+1}]$ is:

$$f_{int}(x) = \sum_{j=0}^{k}l_{i, j, k}(x)f(x_{a_i+j}),$$

where $l_{i,j,k}(x) = \prod_{l=0, \, l\ne j}^k(x-x_{a_i+l})/\prod_{l=0, \, l\ne j}^k(x_{a_i+j}-x_{a_i+l})$ are the Lagrange interpolating polynomials.

Synthesis function.

To derive the spline interpolation basis function from the point of view of $f(x_i)$ I set $f(x_j)=0, \, i\ne j$ and $f(x_i) = 1$. Then I consider the intervals for which $f_{int}$ is nonzero. There are $(k+1)/2$ intervals to the left and $(k+1)/2$ intervals to the right of $x_i$ where it is nonzero. If I enumerate those intervals from left to right, I have $[x_{a_i-1+j}, x_{a_i-1+(j+1)}], \, j = 0, \ldots, k$. Without the re-indexing these are the intervals $[x_{i-(k+1)/2+j}, x_{i-(k+1)/2+(j+1)}]$ whose union spans $[x_{i-(k+1)/2}, x_{i+(k+1)/2}]$. The Lagrange polynomial $l_{a_i-1+j, k-j, k}(x)$ corresponds to the interval $[x_{a_i-1+j}, x_{a_i-1+(j+1)}]$. It should be noted that that as $j$ increases the Lagrange polynomial index $(k-j)$ decreases (so they are flipped).

Ultimately I get the synthesis function:

$$\phi_{k,i}(x) = \begin{cases} l_{a_i-1+j,k-j,k}(x) & x \in [x_{a-1+j},x_{a-1+(j+1)}]\\ 0 & x \in (-\infty, x_{a-1}] \cup [x_{a+k}, \infty) \end{cases}$$

And thus:

$$f_{int}(x) = \sum_i\phi_{k,i}(x)f(x_i).$$

Boundary conditions.

It should be noted that for $0 \leq i < (k+1)/2$ and for $n-(k+1)/2 < i\leq n$ the $\phi_i$ are not well defined since they index points with negative indices, while the available points are $x_0, \ldots, x_n$. One can however extend $x$ appropriately based on the application, for example to be periodic (potentially if the problem is defined on a toroidal domain):

$$x_{n+j} = x_{j-1}, \, x_{-j} = x_{n-j+1}, \, n+1\geq j>0,$$

or mirrored (this is equivalent to requiring that enough derivatives of the interpolant $f_{int}$ are zero at $x_0$ and $x_n$):

$$x_{-j} = x_{j-1}, \, x_{n+j} = x_{n-j+1}, \, n+1\geq j>0,$$

or clamped (this is equivalent to repeating the same point/knot multiple times):

$$x_{-j} = x_0, \, x_{n+j} = x_n, \, j\geq 0.$$

Even degrees.

A similar thing can be done for even degree interpolating basis functions, but by applying the same approach they will be biased either to the left or the right. This is because the even degree Lagrange interpolating polynomial at $[x_i,x_{i+1}]$ needs to pass through an odd number of points $k+1$, so there is no symmetric choice wrt $[x_i,x_{i+1}]$. A workaround would be to consider the average of the interpolating polynomial to the right and to the left as suggested in the answer here Piecewise quadratic interpolation in a symmetric manner. These will still interpolate all of the points that coincide between the two schemes and will only approximate the fringe point on the left and on the right.

Regular grid.

The above ideas can be specialized for a regular grid. On a regular grid the Lagrange interpolating functions $l_{i,j,k}(x)$ will be rescaled and shifted version of $k+1$ canonic ones:

$$l_{j, k}(x) = \frac{\prod_{l=0, l \ne j}^k(x-x_l)}{\prod_{l=0, l \ne j}^k(x_j-x_l)}, \, j = 0, \ldots, k,$$

where the canonic points $x_i$ can be chosen as $x_i = b + ih, i=0, \ldots, k$ for any step $h\ne 0$ and any offset $b\in\mathbb{R}$. Then whenever $l_{i,j,k}(x)$ must be computed whose definition relies on points $x'_i = b' + ih'$, the remapping: $x = b+\frac{h}{h'}(x'-b')$ can be used. Then the synthesis functions $\phi_{k,i}$ can also be written through a single synthesis function $\phi_k$ with a shifted and rescaled argument.

Higher dimensions.

One way to extend the above to higher dimensions is through tensor products. Then Lagrange functions are defined independently for the different grids forming the tensor product.

Continuity.

As mentioned the lagrange synthesis functions do not seem to offer higher continuity unlike B-Splines. On the other hand B-Splines require solving a system of linear equations in order to be interpolating along with achieving $C^{k-1}$. The natural question arises whether continuity lower than $C^{k-2}$ can be achieved with interpolating synthesis functions that require only local information (i.e. do not require a solution of a system of linear equations). At least for $k=3$ the answer seems to be yes, as $C^1$ is achieved with the synthesis function from Keys' paper Cubic convolution interpolation for digital image processing:

$$\phi_3(x) = \begin{cases} \frac{3}{2}|x|^3-\frac{5}{2}x^2+1 & |x| \leq 1 \\ -\frac{1}{2}|x|^3 + \frac{5}{2}x^2 - 4|x| + 2 & 1 \leq |x| \leq 2 \\ 0 & |x|>2 \end{cases}$$

Granted this is derived for a regular grid. It does suggest a way of deriving the synthesis functions with the desired constraints not based on Lagrange or Bernstein bases, but directly with desired constraints on the derivatives at different points and so on.

For a comparison I have added two images (the first is cubic Lagrange and the second is cubic Keys ($C^1$)):

On an unrelated note, funnily the two agree at $\pm 2, \pm \frac{3}{2}, \pm 1, \pm \frac{1}{2}, 0$. The fact that they intersect in the intermediate points makes it so that vertex-centered cubic transfer operators in multigrid could have been derived by sampling either. This is not the case for cell-centered operators however.

Some references.

Some more discussion regarding the Lagrange synthesis functions can be found in Schaum's paper Theory and Design of Local Interpolators (called LF-$(k+1)$ there, see equations (1) and (50)). Finally a whole family of synthesis functions are presented in Muntingh's paper Symbols and exact regularity of symmetric pseudo-splines of any arity (see table 2), however the paper is fairly involved and only the symbols (Fourier transform) are presented, so one has to derive the functions in the spatial domain from that (if not from the generating functions presented there). A similar overview is presented in the paper of Charina et al. Multigrid methods: grid transfer operators and subdivision schemes (see definition 4.1 and example 4.2).