Julien Jerphanion — Personal Website

`tl;dr`

This exposes a way to efficiently sample points from spline curves and 2-sphere homeomorphic splines using linear algebra operations and trivial yet useful caching operations.

First, a definition of spline curves as well as two ways of sampling it is given. Secondly, a definition of 2-sphere homeomorphic splines as well as a matricial scheme to sample several of them efficiently are given.

Finally, some tricks for the implementations are also covered briefly. I think those schemes benefit from a nearly optimal arithmetic intensity.

Notation

In the following, we clarify the following notations:

$[\![n]\!] \triangleq \{0, 1, \cdots, n - 1\}$
As to clarify the manipulation of vectors as rows, for $n \in \mathbb{N}$, $\mathbb{R}^{n}$ is identified to $\mathbb{R}^{1 \times n}$:

$$ \mathbb{R}^n \approxeq \mathbb{R}^{1 \times n} $$

Open or closed spline curves

Expression of the curve $\gamma$

A 2D open or closed Spline is $\gamma$ is defined as follows:

$$ \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]\!],\quad \forall t \in [0,1] \quad \quad \varphi_{\text{per}}(M t - k) = \varphi_{} \circ \text{wrap}_{M}(t,k) $$

Wrapping will be used in all the coming details.

Sampling using matrices

One can efficiently, for $t \in [0, 1]$ 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}^{N-1} \in \mathbb{R}^{N\times M} $$

one gets:

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

Sampling evenly on the curve $\gamma$

A naive 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 2-Sphere homeomorphic Spline 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 the following parameters are considered:

$s\in [0,1]$, for the longitude,
$t\in [0,1]$, for the latitude,

where the following spline coefficients are considered:

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

and where the exponential spline $\phi_a$ for $\alpha \in \left\{\frac{2\pi}{M_t}, \frac{\pi}{M_s - 1}\right\}$ is chosen

\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}

On coefficients and control points

The coefficients shown above contain more than the actual control points.

Indeed, the parametrisation 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):

$$ \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 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:

$$ \left\{ \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} \right. $$

Periodising latitudes using argument wrapping

As in the 2D case, one 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 sums in a sole one, we define $\tilde\Phi(s, t)$ and $\mathbf{\tilde c}$, which are 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]_{l\in [\![-1, M_s]\!], k \in [\![0, M_t ]\!]} \end{array} $$

$$ \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$ and $\varphi$ are fixed? One can proceed as follows:

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 outer 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 with respect to the two last 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 subscripts used for Einstein notations.

About the implementation of $\text{flatten}$

$\text{flatten}$ is implemented using reshaping (e.g np.reshape) which generally is a constant time and memory operation as solely the tensors strides are being modified when reshaping. For more information, see the this comment.

Further work

This Python implementation allows sampling 10 million surface splines’ points in a few second which should be fine for most application.

Still, the slow part of constructing $\mathbf{\Phi}$ can be further speeded up, especially 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?

Finally and more importantly, some time is needed to reorder the implementation, release it and port it to use GPUs.

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

Vectorised implementation of splines

All you need is `numpy.einsum`.