# Vectorised implementation of splines

## All you need is numpy.einsum.

Draft on 25.03.2020
Estimated reading time: ~ 10 min.
Published on 04.10.2021

## 2D open or closed Spline

### Expression of the curve $$\gamma$$

$$\gamma: t \in [0, 1] \longrightarrow \mathbb{R}^2$$
$$\gamma: t \longmapsto \sum_{k=0}^{M-1} \mathbf{c}[k]\ \varphi_{\text{per}}(Mt - k)$$

where $$t$$ is the parameter, where $$\mathbf{c} \in \mathbb{R}^{M \times 2}$$ are the spline coefficients also referred to in this 2D case as control points, and where $$\varphi_{\text{per}}$$ is the periodisation of $$\varphi$$, a given basis function.

### Periodisation using wrapping

Periodisation is needed in the case of closed spline to meet boundary conditions.

One can define a wrapping function $$\text{wrap}_{M}$$ to easily make $$\varphi$$ periodic, that is:

$$\forall k \in [\![M]\!], \forall t \in [0,1] \quad \varphi_{\text{per}}(M t - k) = \varphi_{} \circ \text{wrap}_{M}(t,k)$$

Wrapping is used in this implementation.

### Sampling using matrices

One can efficient sample a point $$\gamma(t) \in \mathbb{R}^2 \approxeq \mathbb{R}^{1\times 2}$$ using matrices multiplication.

One can compute the vector:

$$\phi_M(t) \triangleq \left[\varphi_{\text{per}}(Mt - k)\right]_{k=0}^{M-1} \in \mathbb{R}^M \approxeq \mathbb{R}^{1\times M}$$

And get:

$$\gamma(t) = \phi_M(t)\ \mathbf{c}$$

Thus, for a given number of points $$N$$, one can easily extend this scheme to sample a set $$\mathbf{X}$$ of points:

$$\mathbf X \triangleq [\gamma(t_i)]_{i=0}^{N-1} \in \mathbb{R}^{N\times 2}$$

By computing the matrix:

$$\mathbf{\Phi} \triangleq [\phi_M(t_i)]_{i=0}^{M-1} \mathbb{R}^{N\times M}$$

one gets:

$$\mathbf X = \mathbf{\Phi}\ \mathbf{c}$$

### Sampling evenly on the curve $$\gamma$$

A simple sampling consider a range of $$N$$ evenly spread parameters $$(t_i)_{i=0}^{N-1} \in [0, 1]^N$$, that is:

$$t_i \triangleq \frac{i}{N-1}, \quad \forall i \in [\![N]\!]$$

This generally does not give evenly sampled points $$(\gamma(t_i))_{i=0}^{N-1}$$ on the curve, in the sense that:

$$|| \gamma(t_{i+1}) - \gamma(t_{i})||_2 \approxeq d, \quad \forall i \in [\![N-1]\!]$$

where

$$d\triangleq \frac{L}{N}$$
$$L \triangleq \text{length}_\gamma(1)$$

for:

$$\text{length}_\gamma: t \longmapsto \int_0^t ||\gamma'(t) ||_2 \; \text d t$$

To sample $$(\gamma(t_i))_{i=0}^{N-1}$$ evenly, one can define $$(t_i)_{i=0}^{N-1} \in [0, 1]^N$$ as follows:

$$t_i \triangleq \text{length}_\gamma^{-1}\left(\frac{i}{N-1} L \right), \quad \forall i \in [\![N]\!]$$

In practice, one does not compute $$\text{length}_\gamma^{-1}$$ but binary-searches for preimages of $$\text{length}_\gamma$$ after having it computed on a discretisation of $$[0, 1]$$ using $$N_{\text{oversample}} = c \ N$$ points, for a given integer $$c$$, for instance $$c \triangleq 5$$.

This can be done efficiently for all the points using using numerical integration and a schema with the exact same matricial structure as presented above that considers $$\varphi'$$ instead of $$\varphi$$ for sampling $$\gamma'$$.

## 2-Sphere homeomorphic Spline

### Expression of the surface $$\eta$$

A spline surface in the sphere parametrisation is given by:

$$\eta:[0, 1]^2 \longrightarrow \mathbb{R}^3$$
$$\mathbf{\eta}: (s, t) \longmapsto \sum_{l=-1}^{Ms} \sum_{k=0}^{M_t -1} \mathbf{c}[l, k] \varphi_{\frac{\pi}{M_s - 1}}((M_s - 1)s - l)\ \varphi_{\frac{2\pi}{M_t}, \text{per}}(M_t t - k)$$

where $$s$$ relates to the longitude, $$t$$ to the latitude, where we consider the spline coefficients:

$$\mathbf{c} \in \mathbb{R}^{(M_s + 2)\times M_t \times 3}$$

and exponential spline basis:

\begin{align} \varphi_\alpha:x \longmapsto \begin{cases} \frac{\cos(\alpha|x|)\cos(\frac{\alpha}{2})-\cos(\alpha)}{\left(1-\cos(\alpha)\right)} & 0 \leq |x| < \frac{1}{2}\\ \frac{\left(1-\cos(\alpha(3/2-|x|))\right)}{2\left(1-\cos(\alpha)\right)} & \frac{1}{2} \leq |x| < \frac{3}{2} \\ 0 & \frac{3}{2} \leq |x| \end{cases}. \end{align}

for $$\alpha \in \left\{\frac{2\pi}{M_t}, \frac{\pi}{M_s - 1}\right\}$$.

### On coefficients and control points

The surface is fully determined by $$M_t(M_s-2)+6$$ points:

• the north pole $$\mathbf{c}_\mathrm{N}$$
• the south pole $$\mathbf{c}_\mathrm{S}$$
• the tangents at the north pole $$\mathbf{T}_{1,\mathrm{N}}$$, $$\mathbf{T}_{2,\mathrm{N}}$$
• the tangents at the south pole $$\mathbf{T}_{1,\mathrm{S}}$$, $$\mathbf{T}_{2,\mathrm{S}}$$
• $$(M_s - 2)\times M_t$$ actual control points (excluding $$\mathbf{c_N}$$ and $$\mathbf{c_S}$$, the north and south pole), also referred to as control_points :
$$\mathbf{c}[l, k],\quad (l,k) \in [\![1, M_s - 2 ]\!] \times [\![M_t]\!]$$

Moreover, $$4\times M_t$$ additional control points at the extremities of each open longitude are considered to ensure a correct behaviour at the boundaries are used for continuity and smoothness conditions due to the support size of the exponential spline basis:

$$\mathbf{c}[l, k],\quad (l,k) \in \{-1, 0, M_s - 1, M_s \} \times [\![M_t]\!]$$

Their expression is given, for $$k \in [\![M_t]\!]$$, by:

$$\begin{array} \mathbf{c}[-1,k] &=& \mathbf{c}[1,k] - \frac{\mathbf{T}_{1,\mathrm{N}}c_{M_t}[k] + \mathbf{T}_{2,\mathrm{N}}s_{M_t}[k]}{(M_s-1)\varphi'_{\frac{\pi}{M_s-1}}(1)} \\ \mathbf{c}[0,k] &=& \frac{\mathbf{c}_\mathrm{N} - \varphi_{\frac{\pi}{M_s-1}}(1)(\mathbf{c}[-1,k]+\mathbf{c}[1,k])}{\varphi_{\frac{\pi}{M_s-1}}(0)} \\ \mathbf{c}[M_s,k] &=& \mathbf{c}[M_s-2,k] - \frac{\mathbf{T}_{1,\mathrm{S}}c_{M_t}[k] + \mathbf{T}_{2,\mathrm{S}}s_{M_t}[k]}{(M_s-1)\varphi'_{\frac{\pi}{M_s-1}}(1)}\\ \mathbf{c}[M_s-1,k] &=& \frac{\mathbf{c}_\mathrm{S} - \varphi_{\frac{\pi}{M_s-1}}(1)(\mathbf{c}[M_s-2,k]+\mathbf{c}[M_s,k])}{\varphi_{\frac{\pi}{M_s-1}}(0)} \end{array}$$

### Periodising latitudes using argument wrapping

As in the 2D case, ine can define a wrapping function $$\text{wrap}_{M_t}$$ to make $$\varphi_{\frac{2\pi}{M_t}}$$ periodic, that is

$$\forall k \in [\![M_t]\!], \forall t \in [0,1] \quad \varphi_{\frac{2\pi}{M_t}, \text{per}}(M_t t - k) = \varphi_{\frac{2\pi}{M_t}} \circ \text{wrap}_{M_t}(t,k)$$

Wrapping is used in this implementation.

### Matricial expression of the surface

We can rewrite the expression

$$\mathbf{\eta} (s, t) = \sum_{l=-1}^{Ms} \sum_{k=0}^{M_t -1} \mathbf{c}[l, k] \varphi_{\frac{\pi}{M_s - 1}}((M_s - 1)s - l)\ \varphi_{\frac{2\pi}{M_t}, \text{per}}(M_t t - k)$$

as a matricial product.

To do so, we first define the following matrix:

$$\begin{array} \ \Phi(s, t) &\triangleq& \left[\varphi_{\frac{\pi}{M_s - 1}}((M_s - 1)s - l) \ \varphi_{\frac{2\pi}{M_t}, \text{per}}(M_t t - k)\right]_{(l, k) \in [\![M_s + 2]\!]\times [\![M_t]\!]} &\in& \mathbb{R}^{(M_s + 2)\times M_t} \end{array}$$

This way we have:

$$\mathbf{\eta} (s, t) = \sum_{l=-1}^{Ms} \sum_{k=0}^{M_t -1} \Phi(s, t)[l, k] \ \mathbf{c}[l, k]$$

To combine both sum in one can, we define flattened version of $$\Phi(s, t)$$ and $$\mathbf{c}$$:

$$\begin{array} \ \tilde\Phi(s, t) &\triangleq& \text{flatten}(\Phi(s, t)) &\in & \mathbb{R}^{(M_s + 2)M_t} \approxeq \mathbb{R}^{1 \times (M_s + 2)M_t} \\ &=& \left[\varphi_{\frac{2\pi}{M_t}, \text{per}}(M_t t - k) \varphi_{\frac{\pi}{M_s - 1}}((M_s - 1)s - l)\right]_{i\triangleq (l\ M_t + k) \in [\![(M_s + 2) M_t ]\!]} \end{array}$$

and

$$\begin{array} \quad \mathbf{\tilde c} &\triangleq& \text{flatten}(\mathbf{c}) &\in& \mathbb{R}^{(M_s + 2)M_t \times 3} \\ &=& \left[\mathbf{c}[k,l]\right]_{i\triangleq (l\ M_t + k) \in [\![(M_s + 2) M_t ]\!]} \end{array}$$

So that we get the nice matricial expression:

$$\mathbf{\eta}(s, t) = \tilde\Phi(s, t)\ \mathbf{\tilde c} \in \mathbb{R}^{1\times 3}$$

### Efficiently sampling using the matricial expression

Now, for some sampling parameters $$N_s, N_t$$, how to obtain:

$$\ \mathbf{X} \triangleq \left(\eta(s_i, t_j)\right)_{(i, j) \in {[\![N_s]\!] \times [\![N_t]\!]}} \in \mathbb{R}^{N_sN_t\times 3}$$

?

And how to do it efficiently for several splines when $$N_s, N_t$$ are fixed?

• Compute the tensor resulting from the concatenation of the points $$\Phi$$ matrices:
$$\mathbf{\Phi} \triangleq \left[\Phi(s_i, t_j) \right]_{(i, j) \in {[\![N_s]\!] \times [\![N_t]\!]}} \in \mathbb{R}^{N_s\times N_t\times (M_s + 2)\times M_t}$$
• Flatten it on its two first axis and on its two lasts, that is:
$$\begin{array} \mathbf{\tilde\Phi} &\triangleq& \text{flatten}(\mathbf{\Phi}) \in \mathbb{R}^{N_s N_t\times (M_s + 2)M_t} \\ &=& \left[\tilde\Phi(s_i, t_j) \right]_{p\triangleq(i N_t+ j) \in {[\![N_s N_t]\!]}} \end{array}$$
• Cache $$\mathbf{\tilde\Phi}$$
• Use $$\mathbf{\tilde\Phi}$$ on several splines’ flattened coefficients $$\mathbf{\tilde c}$$:
$$\ \mathbf{X} = \mathbf{\tilde \Phi} \mathbf{\tilde c} \in \mathbb{R}^{N_sN_t \times 3}$$

### Computing $$\Phi$$ efficiently: some tensor wizardry

First note that $$\Phi(s, t) \in \mathbb{R}^{(M_s + 2) \times M_t}$$ can be computed as the outter product of two vectors

$$\Phi(s, t) = \phi_{M_s}^{\text{long}}(s) \otimes \phi_{M_t}^{\text{lat}}(t)$$

where

$$\left\{ \begin{array} \ \phi_{M_s}^{\text{long}}(s) &\triangleq& \left[\varphi_{\frac{2\pi}{M_s -1}}((M_s - 1)s - l)\right]_{l \in [\![M_s + 2]\!]} &\in& \mathbb{R}^{M_s + 2} \\ \phi_{M_t}^{\text{lat}}(t) &\triangleq& \left[\varphi_{\frac{2\pi}{M_t}, \text{per}}(M_t t - k)\right]_{k \in [\![M_t]\!]} &\in& \mathbb{R}^{M_t} \end{array} \right.$$

This can be generalised on all the points $$(i, j) \in {[\![N_s]\!] \times [\![N_t]\!]}$$. Indeed suppose that we have the two following matrices which concatenates the vectors for each point:

$$\left\{ \begin{array} \ \phi_{(N_s, M_s)}^{\text{long}} &\triangleq& \left[\phi_{M_s}^{\text{long}}(s_i)\right]_{i \in [\![N_s]\!]} &\in& \mathbb{R}^{N_s \times (M_s + 2)} \\ \phi_{(N_t, M_t)}^{\text{lat}} &\triangleq& \left[\phi_{M_t}^{\text{lat}}(t_j)\right]_{j \in [\![N_t]\!]} &\in& \mathbb{R}^{N_t \times M_t} \end{array} \right.$$

$$\mathbf{\Phi}$$ can be computed by performing the outer product w.r.t the two axis of those:

$$\mathbf{\Phi} = \left[\phi_{(N_s, M_s)}^{\text{long}}[i,:] \otimes \phi_{(N_t, M_t)}^{\text{lat}}[j,:]\right]_{(i, j) \in {[\![N_s]\!] \times [\![N_t]\!]}} \in \mathbb{R}^{N_s\times N_t\times (M_s + 2)\times M_t}$$

This can be computed in a call of numpy.einsum with the proper subscripts. See this small numpy snippet to understand the implementation.

## About the implementation of $$\text{flatten}$$

$$\text{flatten}$$ is implemented using reshaping (e.g np.reshape) which are constant time and memory operations (solely the tensors strides are changed when reshaping).

## Further work

This Python implementation allows sampling 10 million surface splines’ points in a few second.

I think it’s fine but it might be further speeded-up using parallelization because np.einsum is single-threaded. Might Python transpilers like Numba or Pythran help in this regards as they can translate numpy instructions into efficient machine code?

## References

R. Delgado-Gonzalo, P. Thévenaz, M. Unser, Exponential splines and minimal-support bases for curve representation, Computer Aided Geometric Design, Volume 29, Issue 2, 2012, Pages 109-128, ISSN 0167-8396, https://doi.org/10.1016/j.cagd.2011.10.005