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