# Vectorised implementation of splines

## 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 gets:

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

### Expresion 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 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}
$$

### Perioding 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: the tensor wizardry behind `_get_phi`

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 points each axis:

$$
\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](https://gitlab.ebi.ac.uk/-/snippets/51) 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).

## References