Geometry of Bézier Splines¶
Tip: Launch live version of this tutorial:
It happens frequently that one is interested in the statistical analysis of curves. For example, in order to analyze longitudinal data (e.g., measurements from several subjects at different points in time in a medial study) statistically, hierarchical models (see also this paper) are widely used. By first approximating subject-wise trends (the individual level) and subsequently computing the statistics of these curves (the group level), they account for the fact that intra-subject measurements are highly correlated while inter-subject measurements are not.
In order to analyze Bézier splines thorugh a Riemannian manifold $M$ statistically, we need a Riemannian structure on the space of Bézier curves that go through $M$. Such a structure is implemented in Morphomatics; see A Nonlinear Hierarchical Model for Longitudinal Data on Manifolds for the definitions. In this tutorial we show how it can be used to compute the mean of two Bézier splines in the 2-sphere $S^2$.
%env JAX_PLATFORM_NAME=cpu
import jax
import jax.numpy as jnp
from morphomatics.geom import BezierSpline
from morphomatics.manifold import Bezierfold
from morphomatics.manifold import Sphere
env: JAX_PLATFORM_NAME=cpu
We start by initializing $S^2$ and the manifold of Bézier splines (''Bézierfold'') consisting of splines in $S^2$ with two quadratic segments; we also initialize the Riemannian structure on the Bézierfold.
M = Sphere()
B = Bezierfold(M, 2, 2)
The following lines create two splines as data.
North = jnp.array([0., 0., 1.])
South = jnp.array([0., 0., -1.])
p1 = jnp.array([1., 0., 0.])
o1 = jnp.array([1/jnp.sqrt(2), 1/jnp.sqrt(2), 0.])
om1 = M.connec.exp(o1, jnp.array([0, 0, -.25]))
op1 = M.connec.exp(o1, jnp.array([0, 0, .25]))
q1 = jnp.array([0, 1, 0.])
B1 = BezierSpline(M, [jnp.stack((p1, om1, o1)), jnp.stack((o1, op1, q1))])
z = M.connec.geopoint(o1, North, .5)
p2 = jnp.array([1., 0., 0.])
o2 = M.connec.geopoint(p1, z, .5)
om2 = M.connec.geopoint(p1, z, .4)
op2 = M.connec.geopoint(p1, z, .6)
q2 = z
B2 = BezierSpline(M, [jnp.stack((p2, om2, o2)), jnp.stack((o2, op2, q2))])
Now, we can compute the mean of the two splines and plot the results. Note that the solution is computed iteratively. Therefore, the parameters n, delta, and nsteps can be used to balance speed and accuracy: The larger n and nsteps and the smaller delta, the more accurate (but the slower) the computation will be.
data = jnp.array([B.to_coords(B1), B.to_coords(B2)])
mean = Bezierfold.FunctionalBasedStructure.mean(B, data)[0]
mean = B.from_coords(mean)
Finally, we visualize the two splines (with color gradient) and their mean (solid red).
import numpy as np
import pyvista as pv
# Plotting
# create a tubular surface interpolating points
def tube_mesh(pts, interval=(0,2)):
n = len(pts)-1
curve = pv.PolyData(pts)
curve["time"] = np.linspace(*interval, curve.n_points)
curve.lines = np.c_[np.full(n,2), np.arange(n), np.arange(n)+1]
curve.tube(radius=0.01, inplace=True)
return curve
pl = pv.Plotter(notebook=True)
sphere = pv.Sphere(1)
tube1 = tube_mesh(np.asarray(jax.vmap(B1.eval)(jnp.linspace(0., 2., num=100))))
tube2 = tube_mesh(np.asarray(jax.vmap(B2.eval)(jnp.linspace(0., 2., num=100))))
pl.add_mesh(tube_mesh(np.asarray(jax.vmap(mean.eval)(jnp.linspace(0., 2., num=100)))), color='red')
pl.add_mesh(tube1)
pl.add_mesh(tube2)
pl.add_mesh(sphere)
pl.show_axes()
pl.camera.elevation -= 20
pl.camera.zoom(1.5)
pl.show(jupyter_backend='static')
