We are given a set of data points in the Euclidean space associated with a timestamp $t_i \in \mathbb{Z}$. We would like to compute a composite Bézier curve that interpolates the data points $d_i$ at time $t_i$. In this problem, we try to interpolate four data points $d_i \in \mathbb{S}^2$, the unit sphere. These points are places such that they represent a triangle.
To fit the curve $\mathfrak{B}(t): [0,3] \to \mathbb{S}^2$, you can use our code here, general for manifolds. The possible manifolds are the Euclidean space, the sphere, the special orthogonal group SO(3),... actually all manifolds you can find in Manopt are operational. It is however mandatory to add manopt to your working path. The zip file hereunder contains an old version of Manopt for consistency. You can download the last version of Manopt here (recommended) : www.manopt.org.
Show/hide the code tutorial.
The code
Here is the code for interpolating data with a hybrid composite Bézier curve on the sphere.
cd manopt;
addpath(genpath(pwd));
cd ..;
addpath(genpath([pwd,'/tools']));
phi = [pi/3 ; 4*pi/7 ; 4*pi/7 ; pi/3]  2*pi/12;
theta = [pi/2 ; pi/4 ; 3*pi/4 ; pi/2] + 10*pi/12;
x = cos(theta).*sin(phi); y = sin(theta).*sin(phi); z = cos(phi);
points = [x y z];
data = reshape(points',1,3,4);
n = size(data,3);
t = linspace(0,n1,100);
problem.M = spherefactory(3);
problem.data = data;
problem.control = cp_arnould(problem);
problem = bezier_reconstruction(problem,t);
problem.velocity = bezier_velocity(problem,t);
problem.acceleration = bezier_acceleration(problem,t);
y = squeeze(problem.curve)'; b = squeeze(problem.control)'; p = squeeze(problem.data)';
a = problem.acceleration; v = problem.velocity;
[X,Y,Z] = sphere();
figure;
subplot(2,2,[1 3]); surf(0.95*X,0.95*Y,0.95*Z,...
'FaceAlpha',0.7,'EdgeAlpha',1,...
'FaceColor', [238 197 145]/255);
hold on;
plot3(y(:,1),y(:,2),y(:,3),'b','LineWidth',2);
plot3(b(:,1),b(:,2),b(:,3),'.','Color',[0 0.7 0],'Markersize',30);
plot3(p(:,1),p(:,2),p(:,3),'.r','Markersize',30);
title('Bezier spline');
subplot(222); plot(t(2:end1),v,'b','LineWidth',2); title('velocity');
subplot(224); plot(t(2:end1),a,'b','LineWidth',2); title('acceleration');
Okay, let's dig a bit into that...
Data points
The data points must be stored in a (l,c,n) matrix, where (l,c) is the matrix size of the embedded space. Here, for the sphere, it would be (1,3), as $\mathbb{S}^2$ is embedded in $\mathbb{R}^3$. The points here are draw a triangle on the sphere, and the first and last points coincide.
Hint: You may also call data_points(manifold,type)
to create some predefined test sets. Do not forget to add the path to it (addpath([pwd,'/data_points']);
).
phi = [pi/3 ; 4*pi/7 ; 4*pi/7 ; pi/3]  2*pi/12;
theta = [pi/2 ; pi/4 ; 3*pi/4 ; pi/2] + 10*pi/12;
x = cos(theta).*sin(phi); y = sin(theta).*sin(phi); z = cos(phi);
points = [x y z];
data = reshape(points',1,3,4);
Problem structure preparation
The system is based on a structure called successively by different subfunctions in matlab. Each subfunction will use the elements of differential geometry (the exponential, the logarithm, the distance, etc.) extracted from Manopt and stored in that structure. The names of the different elements of the structure will be very important as they are used by the rest of the code later on.
We prepare first the time parameter t
to be able, later on, to compute the reconstruction spline. Typically, it is a discretization of the domain (here $t \in [0,3]$).
n = size(data,3);
t = linspace(0,n1,100);
Now, we store the crucial information in a structure.
problem.M = spherefactory(3);
problem.data = data;
The structure problem
is supposed to be at the end composed of the following elements (in this specific problem):
problem =
data: [1×3×4 double]
M: [1×1 struct]
control: [1×3×8 double]
curve: [1×3×100 double]
velocity: [98×1 double]
acceleration: [98×1 double]
The data points are stored in the field data
, the control points in control
, the curve in curve
and the manifold in M
. It is important to respect these four names. The other names have less importance.
Bézier curve optimization and reconstruction
The interpolation problem works in two steps:
 find the control points leading the Bézier curve
 reconstruct the actual spline.
Phase 1: Control points generation
The subfunction takes the structure as entry and will return the control points. Store it in the structure for reconstruction.
problem.control = cp_arnould(problem);
There exist different methods to compute the control points:

cp_arnould(problem)
: technique proposed in the paper Arnould et al. 2015. This technique is recommended.

cp_gousenbourger(problem)
: suboptimal technique proposed in the paper Gousenbourger et al. 2014, where the velocity vectors must be specified in the beginning of the method. It works but... well, there is better, so why not using it ? ;).
Phase 2: Reconstruction of the curve
Here also, the subfunction takes the structure as entry and will update it by adding the computed spline. If the structure does not contain the field control
, filled with the control points returned by the control points generation method, the reconstruction is not possible.
problem = bezier_reconstruction(problem,t);
The method is a simple De Casteljau algorithm adapted to compute a hybrid composite Bézier curve.
Optional phase: Acceleration and velocity of the curve
Once the curve is reconstructed, we offer a generic code to evaluate the velocity and the acceleration of the curve. Just use the following piece of code:
a = bezier_velocity(problem,t);
v = bezier_acceleration(problem,t);
Be sure that the field curve
is available in the structure problem
.
The rest is just plotting the curve on a sphere. Try it :) !