Data Fitting on Manifolds with Composite Bézier-Like Curves and Blended Cubic Splines

Abstract

We propose several methods that address the problem of fitting a $C^1$ curve $\gamma$ to time-labeled data points on a manifold. The methods have a parameter, $\lambda$, to adjust the relative importance of the two goals that the curve should meet: being “straight enough” while fitting the data “closely enough.” The methods are designed for ease of use: they only require to compute Riemannian exponentials and logarithms, they represent the curve by means of a number of tangent vectors that grows linearly with the number of data points, and, once the representation is computed, evaluating $\gamma (t)$ at any t requires a small number of exponentials and logarithms (independent of the number of data points). Among the proposed methods, the blended cubic spline technique combines the additional properties of interpolating the data when $\lambda \rightarrow \infty$ and reducing to the well-known cubic smoothing spline when the manifold is Euclidean. The methods are illustrated on synthetic and real data.

Appendix

1.1 Interpolating Bézier Curves: Coefficient Matrices for Control Points Generation

The problem (21) on $\mathcal {M}= \mathbb {R}^m$ is a quadratic function to be optimized with respect to the $n+1$ optimization variables $X = (b_0^+,b_1^-,b_2^-,\dots ,b_n^-)$. The solution of that problem reduces to m independent linear systems $A \cdot X_k = C \cdot D_k$, where $X_k$ is the vector of the $k^\text {th}$ component of the points of X, and $D_k$ is the vector of the $k^\text {th}$ component of the data points $(d_i)_{i=0}^n$. We obtain

$$\begin{aligned}&\int _i^{i+1} \Vert \ddot{\beta }_3(t-i;d_i,b_i^+,b_{i+1}^-,d_{i+1}) \Vert ^2_2 \mathrm {d}t \\&\quad =12 ( -3 d_i b_i^+ + 3 b_i^+ b_i^+ - 3 b_i^+ b_{i+1}^- \\&\qquad + 3 b_{i+1}^- b_{i+1}^- - 3 b_{i+1}^- d_{i+1} + K), \end{aligned}$$

with $i=\lfloor t \rfloor$. K gathers the terms that are independent from the optimization variables. Introducing the differentiability constraints (16) $b_i^+ = 2p_i - b_i^-$, and $\beta _3^i$, the $i^\text {th}$ segment of the composite cubic Bézier curve, one has

$$\begin{aligned} \frac{\partial \beta _3^0}{\partial b_0^+}&= 6 b_{0}^+ - 3 b_{1}^- - 3 d_0 \end{aligned}$$

(42)

$$\begin{aligned} \frac{\partial \beta _3^0}{\partial b_1^-}&= -3 b_{0}^+ + 6 b_{1}^- - 3 d_1 \end{aligned}$$

(43)

and for $i = 1,\dots ,n-1$, $n \ge 2$

$$\begin{aligned} \frac{\partial \beta _3^i}{\partial b_i^-}&= 6 b_{i}^- + 3 b_{i+1}^- - 9 d_i \end{aligned}$$

(44)

$$\begin{aligned} \frac{\partial \beta _3^i}{\partial b_{i+1}^-}&= 3 b_{i}^- + 6 b_{i+1}^- - 6 d_i - 3 d_{i+1}. \end{aligned}$$

(45)

By Definition 2.9, (21) is minimized when these quantities vanish, which yields the linear system $A \cdot X = C \cdot D$, for $n \ge 2$ where (in MATLAB indexing)

$$\begin{aligned}&A(\mathtt{1,1:2}) = \begin{bmatrix} 6&-\,3 \end{bmatrix} \\&A(\mathtt{2,1:3}) = \begin{bmatrix} -\,3&12&3 \end{bmatrix}\\&A(\mathtt{k,k-1:k+1}) = \begin{bmatrix} 3&12&3 \end{bmatrix} \quad k = 3,\dots ,n\\&A(\mathtt{n+1,n:n+1}) = \begin{bmatrix} 3&6 \end{bmatrix}. \end{aligned}$$

and

$$\begin{aligned}&C(\mathtt{1,1}) = 3 \\&C(\mathtt{2,2}) = 12 \\&C(\mathtt{k,k-1:k}) = \begin{bmatrix} 6&12 \end{bmatrix} \quad k = 3,\dots ,n\\&C(\mathtt{n+1,n:n+1}) = \begin{bmatrix} 6&3 \end{bmatrix}. \end{aligned}$$

The third lines in the definition of A and C only hold for $n>2$. All the other entries are equal to zero.

1.2 Fitting Curves : Coefficient Matrices for Control Points Generation

The problem (26) on $\mathcal {M}= \mathbb {R}^m$ is a quadratic function to be optimized with respect to the $2n+2$ optimization variables $X = (p_0,b_0^+,b_1^-,b_1^+,\dots ,b_n^-,p_n)$. As in “Appendix A.1,” the solution of that problem reduces to m independent linear systems $(A_0 + \lambda A_1) \cdot X_k = \lambda C \cdot D_k$. This system depends on the regularization parameter $\lambda > 0$, on X, and on the points $(d_i)_{i=0}^n$ in D.

For $n \ge 4$, the matrices of coefficients $A_0, A_1 \in \mathbb {R}^{(2n+2) \times (2n+2)}$ and $C \in \mathbb {R}^{ (2n+2) \times (n+1)}$ are given by the following sparse matrices.

$A_0$ is given, for $i = 2,\dots ,n-2$, by

$$\begin{aligned}&A_0(\mathtt{1,1:4}) = \begin{bmatrix} 24&-\,36&6&6 \end{bmatrix} \\&A_0(\mathtt{2,1:4}) = \begin{bmatrix} -\,36&72&-\,36&0 \end{bmatrix}\\&A_0(\mathtt{3,1:6}) = \begin{bmatrix} 6&-\,36&48&-\,24&3&3 \end{bmatrix}\\&A_0(\mathtt{4,1:6}) = \begin{bmatrix} 6&0&-\,24&48&-\,33&3\end{bmatrix}\\&A_0(\mathtt{2i+1,2i-1:2i+6}) = \begin{bmatrix} 3&-\,33&48&-\,24&3&3 \end{bmatrix}\\&A_0(\mathtt{2i+2,2i-1:2i+6}) =\begin{bmatrix} 3&3&-\,24&48&-\,33&3 \end{bmatrix}\\&A_0(\mathtt{2n-1,2n-3:2n+2}) =\begin{bmatrix} 3&-\,33&48&-\,24&0&6 \end{bmatrix}\\&A_0(\mathtt{2n,2n-3:2n+2}) = \begin{bmatrix} 3&3&-\,24&48&-\,36&6\end{bmatrix}\\&A_0(\mathtt{2n+1,2n-1:2n+2}) =\begin{bmatrix} 0&-\,36&72&-\,36\end{bmatrix}\\&A_0(\mathtt{2n+2,2n-1:2n+2}) =\begin{bmatrix} 6&6&-\,36&24\end{bmatrix}. \end{aligned}$$

The coefficients of $A_1$ are

$$\begin{aligned}&A_1(\mathtt{1,1}) = \lambda \\&A_1(\mathtt{2,2}) = 0 \\&A_1(\mathtt{2i-1,2i-1:2i}) = \frac{1}{2} \begin{bmatrix} \lambda , \lambda \end{bmatrix} \\&A_1(\mathtt{2n+1,2n+1}) = 0 \\&A_1(\mathtt{2n+2,2n+2}) = \lambda , \end{aligned}$$

for $i = 2,\dots ,n$. Finally, the coefficients of C are given, for $i=2,\dots ,n$, by

$$\begin{aligned}&C(\mathtt{1,1}) = 2\lambda \\&C(\mathtt{2i-1,i}) = \lambda \\&C(\mathtt{2i,i}) = \lambda \\&C(\mathtt{2n+2,n+1}) = 2\lambda . \end{aligned}$$

The other entries are equal to zero.

1.3 Elements of Differential Geometry

Tables 4 and 5 give the explicit formulae used to evaluate the exponential map and the logarithm in this paper. They are implemented in Manopt [7] as a proper factory.

Table 4 Riemannian operators for $\mathbb {S}^{m-1}$ extracted from [31]

Table 5 Riemannian operators for ${\mathrm {SO}}(m)$ extracted from [6]