# Product of Exponentials (POE) Splines on Lie-Groups: Limitations, Extensions, and Application to SO(3) and SE(3)

Andreas Müller

► To cite this version:

Andreas Müller. Product of Exponentials (POE) Splines on Lie-Groups: Limitations, Extensions, and Application to SO(3) and SE(3). International Journal for Numerical Methods in Engineering, 2025, 126 (14), 10.1002/nme.70088. hal-05209022

HAL Id: hal-05209022

https://hal.science/hal-05209022v1

Submitted on 14 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE TYPE

# Product Of Exponentials (POE) Splines on Lie-Groups: Limitations, Extensions, and Application to SO(3) and SE(3)

Andreas Müller

**Correspondence**

Corresponding author Andreas Müller,
Email: a.mueller@jku.at

**Present address**

Johannes Kepler University, Altenberger Str. 69,
4040 Linz

**Abstract**

Existing methods for constructing splines and Bézier curves on a Lie group $G$ involve repeated products of exponentials deduced from local geodesics, w.r.t. a Riemannian metric, or rely on general polynomials. Moreover, each of these local curves is supposed to start at the identity of $G$. Both assumptions may not reflect the actual curve to be interpolated. This paper pursues a different approach to construct splines on $G$. Local curves are expressed as solutions of the Poisson equation on $G$. Therewith, the local interpolations satisfies the boundary conditions while respecting the geometry of $G$. A $k$th-order approximation of the solutions gives rise to a $k$th-order product of exponential (POE) spline. Algorithms for constructing 3rd- and 4th-order splines are derived from closed form expressions for the approximate solutions. Additionally, spline algorithms are introduced that allow prescribing a vector field the curve must follow at the interpolation points. It is shown that the established algorithms, where $k$th-order POE-splines are constructed by concatenating local curves starting at the identity, cannot exactly reconstruct a $k$th-order motion. To tackle this issue, the formulations are extended by allowing for local curves between arbitrary points, rather than curves emanating from the identity. This gives rise to a global $k$th-order spline with arbitrary initial conditions. Several examples are presented, in particular the shape reconstruction of slender rods modeled as geometrically non-linear Cosserat rods.

**KEYWORDS**

Splines, Bézier curves, De Casteljau algorithm, Lie groups, Poisson equation, approximations, product of exponentials, SE(3), SO(3), rigid body motions, Cosserat continua, Magnus expansion

1 | INTRODUCTION

Most systems with practical relevance evolve on Lie groups. Interpolation on manifolds, and Lie groups in particular, are hence relevant in various contexts, such as 3D animation [65, 66], motion planning of robots [29, 40, 61, 68] and unmanned aerial vehicles (UAV) [20, 43], interpolation of motion within numerical integration schemes for multibody systems [5, 27, 59], and the geometric modeling of Cosserat rods [50, 60].

Shoemake [65] was the first who introduced an algorithm for constructing Bézier curves of unit quaternions, i.e. curves on the compact Lie group $Sp(1) \cong SU(2)$, and thus on $SO(3)$, for the purpose of animating spatial rotations. Key element is the spherical linear interpolation (SLERP), which is the representation of a rotation between two given attitudes, $g_0, g_1 \in SO(3)$, by the exponential $\exp(\tau\xi)$ on $SO(3)$, with constant $\xi = \log(g_0^{-1}g_1) \in so(3)$ and path parameter $\tau \in [0, 1]$. The so defined curve corresponds to 1-parameter subgroup of rotations about a constant axis. Moreover, in $SO(3)$, as a compact Lie group, the 1-parameter subgroups are geodesics. This suggests constructing Bézier curves by adopting the classical De Casteljau algorithm on Euclidean vector spaces to compact Lie groups, where geodesics on the group replace straight lines in Euclidean space. A De Casteljau algorithm for constructing Bézier curves on $SO(3)$ was introduced in [55] (this is summarized in appendix B). An algorithm to construct 5th-order Bézier curves for quaternion interpolation, i.e. on $SU(2)$, was presented in [69]. The description of intermediated curves, which are subdivided in course of the De Casteljau algorithm, by an exponential of the form above is valid on any Lie group. This motivated introducing a De Casteljau algorithm on the non-compact Lie group $SE(3)$ of (rigid body)

<footnote n="1">¹Johannes Kepler University, Institute of Robotics, Austria</footnote>
screw motions in [25], where screws are represented as dual quaternions. The intermediate curves (screw motions) described by the exponential are geodesics w.r.t. a left-invariant metric on SE (3). Working in kinematic image space, a cubic interpolation of rigid body motions was introduced in [24]. The Bézier curves on SE (3) were also used to construct ruled surfaces from screw motions [21, 67]. The construction of Bézier curves on symmetric subspaces of SO(3) was addressed in [74]. Clearly, the choice of metric, and its invariance properties, is crucial for constructing meaningful Bézier curves. A thorough mathematical treatment of De Casteljau algorithms on Riemannian spaces can be found in [16], and on Bézier curves and B-splines in [12, 30, 57].

Instead of using the exponential map to describe curves in SE (3) in terms of canonical coordinates, in [64] the Cayley transformation [49] was employed. This gives rise to rational curves. In [58] rational Frenet-Serret motions in SE (3) were introduced by using rational curvature, and the possibility of constructing rotation-minimizing motions was discussed. An algorithm for constructing rational frame motions that minimize rotations was presented in [4]. The groups SO (3) and SE (3) can be modeled as algebraic varieties, e.g. using Study's model [32], and as algebraic groups, which shall provide a framework for developing algorithms to construct rational curves. An exhaustive treatment of rational curves on SO (3) and SE (3) regarded as algebraic groups can be found in [42]. A remarkable result presented in [42] is the generalization of Kempe's Universality theorem, which states that any algebraic curve in the plane can be generated by a mechanical linkage comprising rotational joints only, [41] to homogeneous spaces.

Another approach to the interpolation on Lie groups is adopting the concept of splines defined by piecewise polynomials. In [35, 56], a 2-point interpolation on SO (3) was introduced, where the interpolation between two attitudes is $exp(p(\tau)\xi)$, with a 3rd-order polynomial $p(\tau)$ in $\tau$. Concatenating such 2-point interpolations yields a cubic spline expressed as product of the exponentials of the individual segments. Conceptually the same approach was proposed in [37]. Again, such spline algorithms rely on a distance metric, and an interpolation between two points should be geodesics. A crucial point is the invariance of the particular Riemannian metric when the Lie group is non-compact. This issue was addressed in [17, 18] with the principle of variations. For the special case of interpolation on SE (3), the choice of metric was investigated in [76-79]. Conditions for minimum acceleration and minimum energy curves in SE (3) were derived using a left-invariant metric, e.g. kinetic energy metric. More general discussions on splines in curved spaces can be found in [54], and for general symmetric spaces in [9]. Curves with stationary acceleration on SE (3) were investigated in [63]. It is to be noticed that, due to lack of bi-invariant metric on SE (3), stationary acceleration does not necessarily imply minimum acceleration. Furthermore, the screw motion with canonical coordinates obtained with the logarithm (as above) on SE (3), as used in the screw SLERP [36], is not a geodesic [22]. In conclusion, common to all interpolation methods, is that a certain form of the canonical (exponential) coordinates is deduced from metric properties.

Instead of explicitly relying on metric properties, an alternative approach is to construct interpolation curves as solutions of the Poisson equation on a Lie group $G$, i.e. to reconstruct the curve in the group from a given vector field in the corresponding Lie algebra $g$. The so constructed interpolation naturally accounts for the geometry of $G$, which is encoded in the Poisson equation, but additionally captures the particular motion. In this paper, an approach to interpolation on Lie groups is introduced making use of higher-order solutions of the Poisson equation reported in [47]. It does not rely on geodesics, i.e. on metric properties, as the De Casteljau algorithms in [16, 25, 55, 65], nor does it assume specific polynomial expressions, as the spline algorithms in [35, 56]. A somewhat related publication in this direction is [39], where approximate solutions were used to derive a Newmark integration scheme on SO (3). The approach presented in this paper shares concepts applied to the integration on Lie groups, which can be viewed as interpolation problems when a flow on $G$ is to be reconstructed [46]. Interpolation on SE (3) is an important topic for finite element formulations and Newmark generalized $\alpha$ integration schemes, in general [27, 59]. It should be mentioned that often the interpolation derived from the logarithm on SE (3) is incorrectly regarded as geodesic, e.g. in the finite element formulation in [59].

The general interpolation problem, addressed in this paper, is to find a curve $\bar{h}(t), t \in [0, T]$ in a finite-dimensional connected Lie group $G$ through a given sequence $h_0, h_1, \dots, h_n \in G$, where $\bar{h}(t)$ has required continuity and its derivatives satisfies certain boundary conditions. A more specific problem is the interpolation when additionally samples of a vector field, $v_0, v_1, \dots, v_n \in g$, are given that must be met by the curve. In this paper, spline interpolations are considered that are defined piecewise as

$$ \bar{h}_i(t) = h_0 h_1 \cdots h_{i-1} \exp \xi_i^{[k]}(t), \quad t \in [t_{i-1}, t_i] \qquad (1) $$

where $\bar{h}_i(t) = h_i$, and $\xi_i^{[k]}(t)$ is a curve in the Lie algebra $g$ of $G$ defined as $k$th-order polynomial in $t \in [t_{i-1}, t_i]$. Such a curve will be referred to as a $k$th-order curve in $G$. Since (1) can be written as $\bar{h}_i(t) = h_0 \exp \xi_1^{[k]}(t_1) \cdots \exp \xi_{i-1}^{[k]}(t_{i-1}) \exp \xi_i^{[k]}(t)$, it is referred to as product of exponentials (POE) spline. The special case of cubic splines was reported in [35, 56], where $\xi_i^{[3]}(t)$ are
3rd-order polynomials in $t$, with $\xi_i^{[3]}(t_{i-1}) = 0$. Further, globally defined $k$th-order splines of the form

$$ \bar{h}(t) = \bar{h}_0 \exp \xi^{[k]}(t), \quad t \in [0, T] \qquad (2) $$

are introduced, with non-zero initial value $\xi^{[k]}(0) \neq 0$, and thus $\bar{h}(0) \neq \bar{h}_0$. This enables interpolations of general curves with parameterized initial value. The contributions of this paper are summarized as follows:

* The coordinate functions $\xi_i^{[k]}(t)$ in (1) are derived consistently as $k$th-order approximations of the exact solution of the Poisson equation on $G$ (section 2).

* Closed form expressions for the 3rd- and 4th-order 2-point interpolation with given initial and boundary values on the left trivialized derivative $\bar{v}(t) = \bar{h}(t)^{-1}\bar{h}'(t) \in g$ are presented (section 3). The local order of accuracy of the $k$th-order interpolation is analyzed.

* An algorithm for constructing 3rd- and 4th-order POE-splines through $h_0, h_1, \dots, h_n \in G$ is introduced (section 4.1).

* A 3rd- and 4th-order POE-spline algorithm is presented that admits prescribing values of the vector field $\bar{v}(t) = \bar{h}(t)^{-1}\bar{h}'(t) \in g$ at interpolation points (section 4.2). This enables prescribing terminal conditions on $\bar{h}'$, which is not possible in existing spline formulations.

* A higher-order 2-point interpolation that admits non-zero initial values of the canonical coordinates (section 6).

* An algorithm for constructing globally parameterized $k$th-order splines of the form (2) is proposed (section 7). Such splines can exactly reconstruct a $k$th-order curve in $G$, i.e. a curve $h(t) = h_0 \exp \xi(t)$, with $\xi(t)$ of degree $k$ in $t$. More importantly, they allow for non-zero initial values of the coordinate function $\xi^{[k]}(t)$ in (2). The global spline is thus symmetric in the sense that it admits non-zero terminal as well as initial values.

* The 3rd- and 4th-order splines are applied to the problem of reconstructing the shape of a Cosserat rod described as curve in $SE(3)$, which is particularly important in the area of continuum and soft robotics (examples 1 and 5).

* A De Casteljau algorithm is summarized in appendix B, and results are compared with those of the POE splines.

# 2 | REDUCTION AND RECONSTRUCTION ON LIE GROUPS

## 2.1 | Poisson Equation

Given a curve $g(\tau) \in G$ with (normalized) path parameter $\tau \in [0, 1]$, left translating the vector field $g' \in T_g G$ to the identity $e \in G$ yields the left-invariant vector field $v = g^{-1}g' \in g$. This is known as (kinematic) reduction [8, 45]. Prominent examples in rigid body kinematics and continuum mechanics are rotation group $G = SO(3)$ and the group of rigid body motions $G = SE(3)$. In these cases, $v \in so(3)$ is the body-fixed angular velocity, and $v \in se(3)$ is the body-fixed twist, respectively [53, 62]. The reconstruction problem is to determine $g(t)$ from a given vector field $v(t)$, i.e. to solve the reconstruction equation

$$ g' = gv. \qquad (3) $$

As this equation (in the right-invariant form) was originally derived by Poisson for spatial rotations, it is referred to as (generalized) Poisson equation. It was also reported by Darboux [19], which is why it is also called Poisson-Darboux equation [15].

The solution of (3) can be expressed as

$$ g(\tau) = g_0 \exp \xi(\tau), \quad \tau \in [0, 1] \qquad (4) $$

with $\xi(\tau) \in g$. The so defined $g$ satisfies (3) if and only if $\xi$ satisfies

$$ \xi' = d\exp_{-\xi}^{-1} v \qquad (5) $$

where $d\exp_{\xi}$ is the right-trivialized differential of the exp map [49]. The solution of the ODE system (3) on $G$ is thus reduced to solving the ODE (5) on $g$. Equation (5) is referred to as local reconstruction equation. The solution (5) was reported in the seminal papers [28, 33] and in [70] in terms of a series expansion of the right-trivialized differential. In the following, the series expansion of the inverse differential is used, which was originally reported in [28],

$$ d\exp_{\xi}^{-1} = \sum_{i \ge 0} \frac{B_i}{i!} ad'_{\xi} \qquad (6) $$
where $B_i$ are the Bernoulli numbers. The significance of the (local) reconstruction problem in relation to the interpolation problem lies in the fact that the solution of (3) respectively (5), given specific boundary values, serves as an interpolant. Moreover, since this interpolant satisfies the reconstruction equation, it guarantees consistency with the motion $g(t)$ to be interpolated. This is addressed in the following.

**Remark 1.** The most prominent application of the local reconstruction equations is the family of Munthe-Kaas integration schemes on Lie groups [34, 51, 52], where the original ODE system (3) on $G$ is transcribed to the ODE (5) on $g$, and solved numerically by Runge-Kutta methods (but any other integration can be applied). This is motivated by the fact that (3) must satisfy non-linear invariants whereas (5) must satisfy linear invariants only. The latter are encoded in the vector representation of $g$, and thus automatically preserved during the integration step.

## 2.2 | A kth-Order Approximate Solution of the Poisson Equation

Under the standard assumption $\xi(0) = 0$, it was shown in [47] that the canonical coordinates $\xi(\tau)$ in (4), satisfying (5), admit a series expansion of the form

$$ \xi(\tau) = \sum_{k \ge 0} \frac{\tau^k}{k!} \Phi_k^0 \qquad (7) $$

where the coefficients satisfy the recursive relation

$$ \Phi_k^0 = \mathbf{v}_0^{(k-1)} - (k-1)! \sum_{j=1}^{k-1} \frac{1}{(j-1)!} \mathbf{F}_{kj} \Phi_j^0 \qquad (8) $$

in terms of $\mathbf{v}_0$ and its derivatives, with

$$ \mathbf{F}_{kj} = \sum_{l=1}^{k-j} \frac{(-1)^l}{(l+1)!} \sum_{\substack{i_1, \dots, i_l \in \mathbb{N}_0 \\ i_1 + \dots + i_l + j = k}} \mathrm{ad}_{\mathbf{Y}_{i_l}} \dots \mathrm{ad}_{\mathbf{Y}_{i_l}} \qquad (9) $$

where $\mathbf{Y}_i = \frac{1}{i!} \Phi_i^0$ and $\Phi_0^0 = 0$. The coefficients up to order 5 are

$$ \Phi_1^0 = \mathbf{v}_0, \quad \Phi_2^0 = \mathbf{v}'_0, \quad \Phi_3^0 = \mathbf{v}''_0 + \frac{1}{2} [\mathbf{v}_0, \mathbf{v}'_0], \quad \Phi_4^0 = \mathbf{v}'''_0 + [\mathbf{v}_0, \mathbf{v}''_0] \qquad (10) $$

$$ \Phi_5^0 = \mathbf{v}_0^{(4)} + \frac{3}{2} [\mathbf{v}_0, \mathbf{v}'''_0] + [\mathbf{v}'_0, \mathbf{v}''_0] + \frac{1}{2} [\mathbf{v}'_0, [\mathbf{v}'_0, \mathbf{v}_0]] - \frac{1}{6} [\mathbf{v}_0, [\mathbf{v}_0, \mathbf{v}''_0]] - \frac{1}{6} [\mathbf{v}_0, [\mathbf{v}_0, [\mathbf{v}_0, \mathbf{v}'_0]]] $$

where $[\xi, \eta]$ denotes the Lie bracket of $\xi, \eta \in g$. Throughout the paper, the Lie bracket is also expressed as $[\xi, \eta] = \mathrm{ad}_{\xi}\eta$, where as $\mathrm{ad}: g \times g \to g$ is the adjoint action$^*$ on $g$. Collecting all terms with equal power of $\tau$, the expansion (7) can be written as

$$ \xi(\tau) = \sum_{j=1}^{k} \frac{1}{j!} \tau^j \mathbf{v}_0^{(j-1)} + \sum_{r=3}^{k} \tau^r \mathbf{a}_r \qquad (11) $$

with

$$ \mathbf{a}_r := \sum_{l=2}^{r-1} \sum_{\substack{i_1, \dots, i_l \in \mathbb{N}_0 \\ i_1 + \dots + i_l + l = r}} a_{i_1 \dots i_l} \mathrm{ad}_{\mathbf{v}_0^{(i_1)}} \dots \mathrm{ad}_{\mathbf{v}_0^{(i_{l-1})}} \mathbf{v}_0^{(i_l)}. \qquad (12) $$

The expression (11) is a kth-order extrapolation of the solution of (5) in terms of $\mathbf{v}_0$ and its first $k-1$ derivatives. The non-zero coefficients in (12) are explicitly, for $r = 1, \dots, 5$,

$$ a_3 = \frac{1}{12} [\mathbf{v}_0, \mathbf{v}_0^{(1)}], \quad a_4 = \frac{1}{24} [\mathbf{v}_0, \mathbf{v}_0^{(2)}] \qquad (13) $$

$$ a_5 = \frac{1}{80} [\mathbf{v}_0, \mathbf{v}_0^{(3)}] + \frac{1}{120} [\mathbf{v}_0^{(1)}, \mathbf{v}_0^{(2)}] + \frac{1}{240} [\mathbf{v}_0^{(1)}, [\mathbf{v}_0^{(1)}, \mathbf{v}_0]] - \frac{1}{720} [\mathbf{v}_0, [\mathbf{v}_0, \mathbf{v}_0^{(2)}]] - \frac{1}{720} [\mathbf{v}_0, [\mathbf{v}_0, [\mathbf{v}_0, \mathbf{v}_0^{(1)}]]]. $$

The extrapolation (11) provides a basis for deriving 2-point interpolation formulae of certain approximation orders. It also enables deriving kth-order spline interpolations in a consistent systematic way, in contrast to other approaches, e.g. [35] where general cubic extrapolations, were assumed in the derivation.

<footnote>$^*$ Throughout the paper, the vector representation of Lie algebra elements will be used. For an n-dimensional Lie algebra, $\mathrm{ad}_{\xi}$ is then a $n \times n$ matrix.</footnote>
**Remark 2.** An approximation of the solution of the right Poisson equation, $g' = v_g$, is obtained analogously as above. The solution is expressed as $g(\tau) = \exp \xi(\tau) g_0$, and the right-invariant vector field $v = g'g^{-1}$ satisfies the ODE $\xi' = d\exp_{\xi}^{-1}v$ on $g$. The approximation is obtained with (8), but with the alternating sign $(-1)^l$ in (9) omitted [47]. In case of spatial rotations and rigid body motions, i.e. for Lie groups SO (3) and SE (3), the right-invariant vector field is the spatial representation of angular velocity and twist, respectively.

## 3 | 2-POINT INTERPOLATIONS WITH ZERO INITIAL VALUE $\xi(0) = 0$

**Definition 1.** A kth-order 2-point boundary value interpolation (in canonical coordinates), between $g_0$ and $g_1$ on $G$ is a map

$$ \bar{g}(\tau) = \bar{g}_0 \exp \xi^{[k]}(\tau), \quad \tau \in [0, 1] \qquad (14) $$

satisfying boundary values $\bar{g}(0) = g_0$, $\bar{g}(1) = g_1$, and $\bar{v}^{(s)}(0) = v_0^{(s)}$, $s = 0, \dots, n_0$ and $\bar{v}^{(s)}(1) = v_1^{(s)}$, $s = 0, \dots, n_1$, with $n_0 + n_1 = k - 1$, where $\bar{v} := \bar{g}^{-1}g' \in g$ is the left-trivialized derivative of $\bar{g}$, and $\xi^{[k]}(\tau)$ is a polynomial of degree $k$ in $\tau$.

It is assumed in the following that $\xi^{[k]}(0) = 0$, which implies that $\bar{g}_0 = g(0)$. This is the standard assumption found in the literature. It will be relaxed in section 6 to allow for $\xi^{[k]}(0) \neq 0$. A kth-order $C^{k-1}$ continuos initial value 2-point interpolation between $g_0$ and $g_1$ is derived in the following. The point of departure is the closed form approximate solution (11) of the reconstruction equation (4).

### 3.1 | Initial Value 2-Point Interpolation

It is assumed that the initial values $v_0^{(s)} = v^{(s)}(0)$, $s = 0, \dots, k - 1$ are given, along with the boundary values $g_0$ and $g_1$. At the terminal end, the canonical coordinates are $\bar{\xi} := \log(g_0^{-1}g_1)$. Evaluating the kth-order truncation of the expansion (11) at $\tau = 1$, and solving this for the highest derivative of $v_0$ yields

$$ k! v_0^{(k-1)} = \bar{\xi} - \sum_{j=1}^{k-1} \frac{1}{j!} v_0^{(j-1)} + \sum_{r=3}^{k} a_r. \qquad (15) $$

This is substituted back into (11) giving rise to the kth-order initial value interpolation formula for the canonical coordinates in (14)

$$ \xi^{[k]}(\tau) = \tau^k \bar{\xi} + \sum_{j=1}^{k-1} \frac{1}{j!} (\tau^j - \tau^k) v_0^{(j-1)} + \sum_{r=3}^{k-1} (\tau^r - \tau^k) a_r \qquad (16) $$

with $a_r$ in (12). In particular, the 3rd-, 4th-, and 5th-order initial value interpolation formulae are

$$ \xi^{[3]}(\tau) = \tau^3 \bar{\xi} + (\tau - \tau^3) v_0 + \frac{1}{2} (\tau^2 - \tau^3) v'_0 \qquad (17) $$

$$ \xi^{[4]}(\tau) = \tau^4 \bar{\xi} + (\tau - \tau^4) v_0 + \frac{1}{2} (\tau^2 - \tau^4) v'_0 + \frac{1}{6} (\tau^3 - \tau^4) (v_0'' + \frac{1}{2} [v_0, v_0']) \qquad (18) $$

$$ \xi^{[5]}(\tau) = \tau^5 \bar{\xi} + (\tau - \tau^5) v_0 + \frac{1}{2} (\tau^2 - \tau^5) v'_0 + \frac{1}{6} (\tau^3 - \tau^5) (v_0'' + \frac{1}{2} [v_0, v_0'']) + \frac{1}{12} (\tau^4 - \tau^5) (v_0''' + [v_0, v_0'']) \qquad (19) $$

The boundary conditions $\xi^{[k]}(0) = 0$, $\xi^{[k]}(1) = \bar{\xi}$, and the initial condition on $\bar{v}_0 = \bar{g}^{-1}g'_0|_{\tau=0}$ and on its $k-1$ derivatives are satisfied by construction.

**Remark 3.** It was shown in [56] that, on SO(3), that (17) is the minimum energy curve if $v_0$ and $v'_0$ are zero. If the initial velocity and its derivatives are zero, all kth-order interpolations describe the same path (a 1-parameter subgroup of $G$) but the trajectories differ (kth-order polynomials in $\tau$).

**Lemma 1.** The kth-order initial value 2-point interpolation formula (16) has local error

$$ \tau^k \left( \bar{\xi} - \sum_{j=1}^{k} \frac{1}{j!} v_0^{(j-1)} - \sum_{r=3}^{k} a_r \right) + O(\tau^{k+1}). \qquad (20) $$
Proof. Let $g(\tau) = g_0 \exp \xi(\tau)$ be the curve in $G$ to be interpolated, with $\xi(0) = 0$, and $g^{[3]}(\tau) = \bar{g}_0 \exp \xi^{[k]}(\tau)$ be the interpolant, with $\bar{g}_0 = g_0$. The local error is introduced as $\log(g^{-1}g^{[3]}_0)$. The expansion (11) of $g(\tau)$ with $k \to \infty$ and the $k$th-order interpolation (16), the local error becomes

$$ \begin{aligned} \log(g^{-1}g^{[3]}) &= \xi^{[k]}(\tau) - \xi(\tau) + \frac{1}{2}[\xi^{[k]}(\tau), \xi(\tau)] + \dots \\ &= \tau^k \left( \bar{\xi} - \sum_{j=1}^k \frac{1}{j!} \mathbf{v}_0^{(j-1)} - \sum_{r=3}^k \mathbf{a}_r \right) - \sum_{j>k} \frac{1}{j!} \tau^j \mathbf{v}_0^{(j-1)} - \sum_{r>k} \tau^r \mathbf{a}_r \\ &\quad + \frac{1}{2} \tau^k \left[ \bar{\xi} - \sum_{j=1}^{k-1} \frac{1}{j!} \mathbf{v}_0^{(j-1)} - \sum_{r=3}^{k-1} \mathbf{a}_r, \xi(\tau) \right] + \frac{1}{2} \left[ \sum_{j=1}^{k-1} \frac{1}{j!} \tau^j \mathbf{v}_0^{(j-1)} + \sum_{r=3}^{k-1} \tau^r \mathbf{a}_r, \sum_{j \ge k} \frac{1}{j!} \tau^j \mathbf{v}_0^{(j-1)} + \sum_{r \ge k} \tau^r \mathbf{a}_r \right] \\ &\quad + \frac{1}{2} \left[ \sum_{j=1}^{k-1} \frac{1}{j!} \tau^j \mathbf{v}_0^{(j-1)} + \sum_{r=3}^{k-1} \tau^r \mathbf{a}_r, \sum_{j=1}^{k-1} \frac{1}{j!} \tau^j \mathbf{v}_0^{(j-1)} + \sum_{r=3}^{k-1} \tau^r \mathbf{a}_r \right] \\ &= \tau^k \left( \bar{\xi} - \sum_{j=1}^k \frac{1}{j!} \mathbf{v}_0^{(j-1)} - \sum_{r=3}^k \mathbf{a}_r \right) + O(\tau^{k+1}). \end{aligned} $$

The Lie bracket term in the third line vanishes, and the remaining terms are of order higher than $k$ in $\tau$ leading to the local error (20). $\square$

**Lemma 2.** The $k$th-order initial value 2-point interpolation exactly reconstructs a $k$th-order curve defined on a 1-parameter subgroup of $G$.

Proof. A $k$th-order 1-parameter subgroup is defined by $g(\tau) = g_0 \exp(p(\tau)\xi)$, with constant $\xi \in g$ and a polynomial $p(\tau) = c_0 + c_1\tau + \dots + c_k\tau^k$ of degree $k$. The velocity and derivatives are thus $\mathbf{v}_0^{(s)} = p^{(s+1)}\xi$. The polynomial satisfies $p(0) = 0$ and $p(1) = 1$, thus $c_0 = 0$ and $c_1 + \dots + c_k = 1$. The derivatives of $p$ at $\tau = 0$ are $p^{(s)}(0) = s!c_s$. The 2-point interpolation is then

$$ \xi^{[k]}(\tau) = \tau^k \bar{\xi} + \sum_{j=1}^{k-1} \frac{1}{j!} (\tau^j - \tau^k) \mathbf{v}_0^{(j-1)} + \sum_{r=3}^{k-1} (\tau^r - \tau^k) \mathbf{a}_r = \tau^k \bar{\xi} + \sum_{j=1}^{k-1} c_j (\tau^j - \tau^k) \xi = (c_1 \tau + c_2 \tau^2 + \dots + c_k \tau^k) \xi = p(\tau) \xi $$

hence exactly reproduce the $k$th-order curve. $\square$

The following is a direct consequence of the De Casteljau algorithm 9 in appendix B for constructing Bézier curves on Lie groups. The proof follows immediately from the fact the log and exp map cancel (within their uniqueness domains), and that the Lie group algorithm does reduces to the classical De Casteljau algorithm on a vector space, which can exactly reconstruct a $k$th-order curve with correctly chosen control points.

**Proposition 1.** A Bézier curve of order $k$ constructed with the De Casteljau algorithm 9 (appendix B) on $G$ exactly reconstructs a $k$th-order curve defined on a 1-parameter subgroup of $G$, if the control points (in terms of $\tau$) are the Bernstein basis.

## 3.2 Boundary Value 2-Point Interpolation

Now it is assumed that, in addition to some initial values, the terminal value $\mathbf{v}_1$ is given. The general case, when derivatives of $\mathbf{v}_1$ are prescribed, is not addressed here.

A 2-point interpolation is derived from the initial value interpolation (16). To this end, the derivative of the $k$th-order interpolation (16) is evaluated at $\tau = 1$, which yields, after some manipulation,

$$ \begin{aligned} \frac{d^s}{d\tau^s} \xi^{[k]} (1) = \frac{k!}{(k-s)!} \bar{\xi} - \frac{s}{(k-s)!} \mathbf{v}_0^{(k-2)} + \sum_{j=s}^{k-2} \frac{1}{(j-s)!} \mathbf{v}_0^{(j-1)} - \sum_{j=1}^{k-2} \frac{k!}{j!(k-s)!} \mathbf{v}_0^{(j-1)} \\ + \sum_{r=\max(3,s)}^{k} \frac{r!}{(r-s)!} \mathbf{a}_r - \sum_{r=3}^{k} \frac{k!}{(k-s)!} \mathbf{a}_r \end{aligned} \quad (21) $$
for $1 \le s \le k - 1$. The derivative of $\xi(\tau)$ is also obtained exactly from (5) as

$$ \frac{d^s}{d\tau^s} \xi(\tau) = \Phi_s(\xi, \mathbf{v}, \mathbf{v}', \dots, \mathbf{v}^{(s-1)}) \qquad (22) $$

where

$$ \Phi_k(\xi(\tau), \mathbf{v}(\tau), \mathbf{v}'(\tau), \dots, \mathbf{v}^{(s-1)}(\tau)) := \sum_{i=1}^{k} \Phi_{ki}(\xi(\tau)) \mathbf{v}^{(i-1)}(\tau), \quad k \ge 1 \qquad (23) $$

$$ \text{with } \Phi_{ki}(\xi(\tau)) := \binom{k-1}{i-1} \frac{d^{k-i}}{d\tau^{k-i}} \mathrm{dexp}_{-\xi(\tau)}^{-1} \qquad (24) $$

Evaluated at $\tau = 1$ yields $\frac{d^s}{d\tau^s} \xi(1) = \Phi_s(\bar{\xi}, \mathbf{v}_1, \mathbf{v}'_1, \dots, \mathbf{v}^{(s-1)}_1)$. The latter depends on $\mathbf{v}_1$ and its derivatives, while (21) depends on $\mathbf{v}_0$ and its derivatives. Equating the two at $\tau = 1$, and solving (21) for the highest derivative $\mathbf{v}_0^{(k-2)}$ yields

$$ \mathbf{v}_0^{(k-2)} = \frac{k!}{s} \bar{\xi} - \frac{(k-s)!}{s} \Phi_s(\bar{\xi}, \mathbf{v}_1, \mathbf{v}'_1, \dots, \mathbf{v}^{(s-1)}_1) + \sum_{j=s}^{k-2} \frac{(k-s)!}{(j-s)!s} \mathbf{v}_0^{(j-1)} - \sum_{j=1}^{k-2} \frac{k!}{j!s} \mathbf{v}_0^{(j-1)} + \sum_{r=\max(3,s)}^{k} \frac{r!(k-s)!}{(r-s)!s} \mathbf{a}_r - \sum_{r=3}^{k} \frac{k!}{s} \mathbf{a}_r. \qquad (25) $$

Thus $\mathbf{v}_0^{(k-2)}$ can be expressed in terms of $\mathbf{v}_0, \dots, \mathbf{v}_0^{(k-3)}, \mathbf{v}_1, \dots, \mathbf{v}_1^{(s-1)}$. In the following, only the case $s = 1$ is considered, in which the highest derivative $\mathbf{v}_0^{(k-2)}$ is replaced by the terminal value $\mathbf{v}_1$. This seems to be most relevant for applications.

For the 3rd-order interpolation, i.e. $k = 3, s = 1$, relation (25) yields

$$ \mathbf{v}'_0 = 6\bar{\xi} - 2\mathrm{dexp}_{-\xi}^{-1}\mathbf{v}_1 - 4\mathbf{v}_0. \qquad (26) $$

Inserting this result back into (17) yields the **3rd-order boundary value interpolation**

$$ \xi^{[3]}(\tau) = (3\tau^2 - 2\tau^3)\bar{\xi} + \tau(\tau - 1)^2 \mathbf{v}_0 + (\tau^3 - \tau^2) \mathrm{dexp}_{-\xi}\mathbf{v}_1 \qquad (27) $$

with prescribed initial and terminal velocity. Relation (25) evaluated for the 4th-order interpolation ($k = 4, s = 1$) yields

$$ \mathbf{v}'_0 = 24\bar{\xi} - 6\mathrm{dexp}_{-\xi}^{-1}\mathbf{v}_1 - 18\mathbf{v}_0 - 6\mathbf{v}'_0 - \frac{1}{2}[\mathbf{v}_0, \mathbf{v}'_0] \qquad (28) $$

and, along with (18), leads to the the **4th-order boundary value interpolation**

$$ \xi^{[4]}(\tau) = (4\tau^3 - 3\tau^4)\bar{\xi} + \tau(\tau - \tau^2)(1 + 2\tau)\mathbf{v}_0 + (\tau^4 - \tau^3)\mathrm{dexp}_{-\xi}\mathbf{v}_1 + \frac{1}{2}\tau^2(1 - \tau)^2\mathbf{v}'_0 \qquad (29) $$

with prescribed initial velocity and derivative, and terminal velocity. In the same way the **5th-order boundary value interpolation** ($k = 5, s = 1$), with additionally prescribed initial value $\mathbf{v}'_0''$, is obtained as

$$ \xi^{[5]}(\tau) = (5\tau^4 - 4\tau^5)\bar{\xi} + (\tau - 4\tau^4 + 3\tau^5)\mathbf{v}_0 + \frac{1}{2}(\tau^2 - 3\tau^4 + 2\tau^5)\mathbf{v}'_0 + \frac{1}{6}(\tau^3 - 2\tau^4 + \tau^5)(\mathbf{v}_0'' + \frac{1}{2}[\mathbf{v}_0, \mathbf{v}'_0]) + (\tau^5 - \tau^4)\mathrm{dexp}_{-\xi}\mathbf{v}_1 \qquad (30) $$

**Remark 4.** It was shown in [79] that the 3rd-order approximation (27) is exactly the minimum acceleration curve when interpolating rigid body motions on the Lie group SE(3). This was also observed in [56] for the special case of zero initial and terminal velocities, $\mathbf{v}_0 = \mathbf{v}_1 = \mathbf{0}$. In the latter case, the minimum acceleration curve (27) reduces to the 1-parameter motion $\xi^{[3]}(\tau) = (3\tau^2 - 2\tau^3)\bar{\xi}$, which is a geodesic w.r.t. the left-invariant metric on SE(3).

**Example 1.** The generation of curves in the Lie group SE(3) = SO(3) × $\mathbb{R}^3$ of isometric and orientation preserving transformations in Euclidean 3-space is of particular practical relevance. Then, $\xi \in \mathbb{R}^6 \cong \mathrm{se}(3)$ is the screw coordinate vector that generates a frame transformation $g \in \mathrm{SE}(3)$ via the exp map (A6). Besides trajectory planning of UAVs [20, 43], the shape reconstruction of flexible rods, modeled as Cosserat continua, is an application of particular importance for soft robots [7, 11, 31, 50]. The rod's cross section performs spatial motions parameterized by the normalized arc length $\tau \in [0, 1]$. The spatial pose of the rod's cross section is represented by a rigid frame transformation $h(\tau) \in \mathrm{SE}(3)$. A Cosserat rod is thus a one-dimensional continuum, represented as a curve in SE(3). The deformation field in body-fixed representation is then defined by the Poisson equation (3) as $\mathbf{v} = g^{-1}\mathbf{g}'$ [11]. As an example, the following problem is addressed: Given the pose $g_0 = g(0)$ and $g_1 = g(1)$ at the start and
terminal end, and the deformation/strain $v_0 = v(0)$ and $v_1 = v(1)$ at the both ends, compute the shape of the rod in a static equilibrium, i.e. the curve $g(t)$ in $SE(3)$ for which the elastic potential is stationary. This is a typical problem for dual arm robotic manipulation of deformable slender objects [1, 2], and for continuum robots [71, 75]. The static equilibrium is governed by the Euler-Poincaré equation [11]

$$ \Lambda' - ad_{\mathbf{v}}^{T} \Lambda = \mathbf{W} \qquad (31) $$

where $\Lambda(\tau) \in se^*$ (3) is the stress, and $\mathbf{W}(\tau) \in se^*$ (3) represents distributed loads along the beam. Assuming linear elastic constitutive relations, the stress is related to the strain as $\Lambda = K(v - v^0)$ with the stiffness matrix $K$. Therein, $v^0$ is in correspondence to the geometry of the undeformed beam. For an initially straight beam, where the x-axis of the cross sectional frame is aligned with the beam center line, it is $v^0 = [0, 0, 0, 1, 0, 0]^T$. The cross section frame is located at the center of the cross section and aligned with the principal axes. Then the stiffness matrix is $K(\tau) = \text{diag}(GJ_x, EI_{yy}, EI_{zz}, EA, GA, GA)$, with Young's modulus $E$, shear stiffness $G$, the cross section area $A$, and the second area moments $I_{xx}, I_{zz}$, and polar moment $J_x$. Inserting the constitutive relation into the Euler-Poincaré equation (31) yields the ODE on $se(3)$ for the deformation field $v$

$$ v' = v_0' + K^{-1}(ad_{\mathbf{v}}^{T}K - K')(v - v^0) + K^{-1}\mathbf{W}. \qquad (32) $$

This ODE can be solved with given initial/terminal value of $v$ along with the kinematic equation (3). The boundary deformation is obtained from (32) as $v_0 = v^0 - K^{-1}(0)\mathbf{W}(0)$ and $v(1) = v^0 - K^{-1}(1)\mathbf{W}(1)$, where $\mathbf{W}(0)$ and $\mathbf{W}(1)$ are concentrated loads at the start and terminal end, respectively. Notice that the above equations are normalized.

The static equation (31), respectively (32), is only used in this example to provide a reference solutions $\xi(\tau)$, $v(\tau)$. These solutions are compared with the 3rd- and 4th-order interpolations obtained with $v_0$ and $v_1$. The exponential map on $SE(3)$ and its right-trivialized derivative expd in closed form, and the definition of Lie bracket on $se(3)$ are summarized in appendix A. For computations, a clamped slender rubber beam of 100 mm length, a constant $8 \times 8 \text{ mm}^2$ square cross section, and material parameters $E = 10 \text{ MPa}$, $G = 0.3 \text{ MPa}$ is considered. Distributed loads $\mathbf{W}$ are neglected, thus only the wrench $\mathbf{W}(0)$ at the start, and $\mathbf{W}(1)$ at the terminal end due to the geometric fixation are present. The ODE (32) and (3) are solved with the Munthe-Kaas method using a RK4 integration scheme [26, 51] to obtain a reference solutions $\xi(\tau)$ and $v(\tau)$. The terminal wrench at $\tau = 0$ is prescribed as $\mathbf{W}(0) = [0, -0.0293, -0.1277, 0.0977, 0.0665, -0.1950]$. From the latter, the deformation at the start end is obtained as $v_0 = [0, 0.85797, 3.71876, 0.99999, -0.00035, 0.00102]$. With $\xi(0) = 0$ and $v(0) = v_0$, the ODE system is solved as initial value problem, which yields the terminal pose of the cross section, relative to its initial pose (Fig. 2),

$$ g_1 = \begin{bmatrix} \mathbf{R}_1 & \mathbf{r}_1 \\ \mathbf{0} & 1 \end{bmatrix} \in SE(3), \quad \mathbf{R}_1 = \begin{bmatrix} -0.51519 & -0.42868 & -0.74217 \\ -0.69850 & -0.29163 & 0.65339 \\ -0.49653 & 0.85509 & -0.14923 \end{bmatrix} \in SO(3), \quad \mathbf{r}_1 = \begin{bmatrix} -0.14518 \\ 0.29880 \\ -0.42338 \end{bmatrix} \in \mathbb{R}^3 $$

and $v_1 = [2.92969, -0.50391, 1, 0.00119, 0]$ and $v_0' = [0, 5.71458, 1.94894, -0.00006, 0.00189, -0.00044]$.

The 3rd- and 4th-order approximation are computed with the 2-point interpolations (27) and (29), respectively. The deviation from the reference solution is computed as $\varepsilon^{[k]}(\tau) = \log(\exp(-\xi(\tau))\exp(\xi^{[k]}(\tau)))$, $k = 3, 4$, where $\varepsilon^{[k]} = [x^{[k]}, y^{[k]}] \in \mathbb{R}^6$ is the screw coordinate vector. As there is no bi-invariant metric on $SE(3)$, the deviation is quantified separately for the $SO(3)$ and $\mathbb{R}^3$ part, as $\varepsilon_r^{[k]}(\tau) = \|x^{[k]}(\tau)\|$ and $\varepsilon_p^{[k]}(\tau) = \|y^{[k]}(\tau)\|$, respectively. Fig. 1a) shows the results obtained with (27) and (29). For completeness, the actual shape is shown in Fig. 2.

Finally, the length of the segment over which the shape is interpolated is varied to show the order of interpolation error. To this end, the terminal point of the interpolation is located at $0 < T \le 1$, i.e. the segment length is $T$, and the terminal values are computed as $\bar{\xi} = \xi(T)$ and $v_1 = v(T)$, $v_1' = v'(T)$. The error measures $\varepsilon_r$ and $\varepsilon_p$ are evaluated at the center of the segment. Fig. 1b) shows values of $\varepsilon_r(T/2)$ and $\varepsilon_p(T/2)$ for $T = 0.01, \dots, 1$. The error decreases with order 3 and 4, respectively, with the segment length in accordance with Lemma 1. The result when interpolation the displacement field of a rod with rectangular cross section is shown in Fig. 3. The error is similar to that of the rod with square cross section in Fig. 1.

## 4 | PIECEWISE PARAMETERIZED POE-SPLINES

A kth-order spline interpolating through $h_0, h_1, \dots, h_n \in G$ is constructed by concatenating 2-point initial value interpolations, defining n spline segments. The spline curve is described by the global path parameter $t \in [0, T]$. In segment i, a kth-order interpolation (16) between $h_{i-1}$ and $h_i$ is introduced. The range of the path parameter in segment i is $t \in [t_{i-1}, t_i]$. A normalized
local path parameter for segment i is introduced as $\tau_i(t) = (t - t_{i-1}) / (t_i - t_{i-1}) \in [0, 1]$. The local coordinates in segment i are denoted with $\xi_i^{[k]}(\tau_i)$. It is assumed that the local coordinates are zero at the begin of a segment, i.e. $\xi_i^{[k]}(0) = 0$.

**Definition 2.** A kth-order C^p-continuos POE-spline interpolating through $h_0, h_1, \dots, h_n \in G$ is a map defined recursively in segment i as

$$ \bar{h}_i(t) = h_{i-1} \exp \xi_i^{[k]}(\tau_i(t)), \quad t \in [t_{i-1}, t_i] \qquad (33) $$

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** a) Error when interpolating a curve in SE(3), representing a flexible rod, with given initial and terminal strain, using the 3rd- and 4th-order interpolation. b) Interpolation error at the midle of the flexible rod as function of segment length T.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**FIGURE 2** Shape of the flexible rod as result of a) the 3rd- and b) the 4th-order interpolation. The exact numerical solution of the static balance equation (32) is shown in gray color as a reference (partially covered by the approximate solutions).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**FIGURE 3** Shape reconstruction of a flexible rod with rectangular cross section. a) 3rd- order (red) and b) 4th-order interpolation (green) are superposed to the numerically exact solution (gray).</figcaption></figure>
satisfying the matching condition $\bar{h}_i(t_i) = h_i$, the continuity conditions $\bar{v}_i^{(s)}(0) = \bar{v}_{i-1}^{(s)}(1)$, $s = 0, \dots, p-1$, and the initial conditions $\bar{v}_1^{(s)}(0) = v_0^{(s)}$, $s = 0, \dots, p-1$ on the vector fields $\bar{v}_i(\tau_i(t)) := \bar{h}_i(t)^{-1} \bar{h}_i'(t)$ and its derivatives $((\cdot)')$ denotes derivative w.r.t. $t$.

**Remark 5.** The map (33) can be written as

$$ \bar{h}_i(t) = h_0 \exp \bar{\xi}_1 \dots \exp \bar{\xi}_{i-1} \exp \xi_i^{[k]}(\tau_i(t)), \quad t \in [t_{i-1}, t_i] \qquad (34) $$

with $\bar{\xi}_i = \log(\bar{h}_{i-1}^{-1} \bar{h}_i)$. This is why (33) is referred to as POE-spline. A POE-spline relies on a *local interpolation* on G around the identity with local coordinates $\xi_i^{[k]}$.

**Remark 6.** In all POE-spline formulations reported in the literature [35,56], it is presumed that $\xi_i^{[k]}(0) = 0$. That is, $\xi_i^{[k]}$ are local coordinates on G used to reconstruct the 'incremental' motion within a segment, starting at the identity of G. This concept is also employed in Munthe-Kaas integration schemes on Lie groups [34, 51, 52]. The consequences of the fact that the spline is constructed piecewise in terms of local coordinates for the n segments separately will be discussed in section 5.

## 4.1 | A kth-Order $C^{k-1}$ Continuos POE-Spline

The matching conditions are equivalent to $\exp \xi_i^{[k]}(1) = \bar{h}_{i-1}^{-1} \bar{h}_i$, which are satisfied by construction of the 2-point interpolation. It remains to impose the continuity conditions. The 2-point interpolation in segment i needs initial values $\bar{v}_i^{(s)}(0)$, $s = 1, \dots, k-1$ (replacing $v_0^{(s)}$ in (16)) that match those at the end of segment i. The vector field obtained from the interpolation in segment i-1 is

$$ \bar{v}_{i-1}(\tau_{i-1}) = \bar{h}_{i-1}(t)^{-1} \bar{h}_{i-1}'(t) = \frac{1}{t_i - t_{i-1}} \exp_{-\xi_{i-1}^{[k]}(\tau_{i-1})} (\xi_{i-1}^{[k]'}(\tau)) \qquad (35) $$

noting that $\frac{d}{dt} = \frac{1}{t_i - t_{i-1}} \frac{d}{d\tau}$. At the end of segment i-1 this is $\bar{v}_{i-1}(1) = \frac{1}{t_i - t_{i-1}} \exp_{-\bar{\xi}_i}(\xi_{i-1}^{[k]'}(1))$, which delivers the initial value for the 2-point interpolation in segment i. The higher-order continuity conditions are

$$ \bar{v}_i^{(s)}(0) = \frac{1}{(t_i - t_{i-1})^{s+1}} \frac{d^s}{d\tau^s} \exp_{-\xi_{i-1}^{[k]}(\tau)} (\xi_{i-1}^{[k]'}(\tau)) \Big|_{\tau=1} \qquad (36) $$

with $\xi_{i-1}^{[k]}(1) = \bar{\xi}_{i-1}$.

For a 3rd-order spline with 2nd-order continuity, these relations are

$$ \bar{v}_i(0) = \frac{1}{t_i - t_{i-1}} \exp_{-\bar{\xi}_{i-1}} (\xi_{i-1}^{[k]'}(1)) \qquad (37) $$

$$ \bar{v}_i'(0) = \frac{1}{(t_i - t_{i-1})^2} \left( \exp_{-\bar{\xi}_{i-1}} (\xi_{i-1}^{[k]''}(1)) - (\mathrm{D}_{-\bar{\xi}_{i-1}} \exp)(\xi_{i-1}^{[k]'}(1))(\xi_{i-1}^{[k]'}(1)) \right). \qquad (38) $$

A 4th-order spline with 3rd-order continuity, must satisfy the additional condition

$$ \bar{v}_i''(0) = \frac{1}{(t_i - t_{i-1})^2} \left( \exp_{-\bar{\xi}_{i-1}} (\xi_{i-1}^{[k]''}(1)) - 2(\mathrm{D}_{-\bar{\xi}_{i-1}} \exp)(\xi_{i-1}^{[k]'}(1))(\xi_{i-1}^{[k]''}(1)) \right. \\ \left. -(\mathrm{D}_{-\bar{\xi}_{i-1}} \exp)(\xi_{i-1}^{[k]''}(1))(\xi_{i-1}^{[k]'}(1)) + (\mathrm{D}_{-\bar{\xi}_{i-1}}^2 \exp)(\xi_{i-1}^{[k]'}(1))(\xi_{i-1}^{[k]'}(1))(\xi_{i-1}^{[k]'}(1)) \right). \qquad (39) $$

The derivatives of $\xi_{i-1}^{[k]}(\tau_i)$ are obtained from (16), particularly from (17) and (18) for $k = 3, 4$, respectively. With these continuity conditions, it is straightforward to derive the 3rd-order $C^2$ continuos spline Algorithm 1, and the 4th-order $C^3$ continuos spline Algorithm 2 below. An equivalent 3rd-order algorithm was presented in [35,56] for the special case of POE-splines on SO(3), where 2-point interpolations are derived as minimum energy curves. That is, the spline algorithm in [35,56], constructed by combining minimum energy curves, is recovered as the POE-spline constructed from the 2-point interpolations derived from 3rd-order approximate solutions of the Poisson equation. The computation effort of the POE spline algorithm grows linear with the number of points, i.e. it has linear complexity.
Algorithm 1: 3rd-Order C² continuos POE-Spline

1. Input:

• Knot points $t_0, t_1, \dots, t_N$

• Samples $h_0, h_1, \dots, h_n$

• Initial values $v_0, \dot{v}_0$

3. Computation of Spline Coefficients:

FOR $i = 1, \dots, N$ DO

• $T_i := t_i - t_{i-1}$

• $\delta_i := T_i/T_{i-1}$

• $\bar{\xi}_i := \log(h_{i-1}^{-1}h_i)$

• $\alpha_{i-1} := \delta_i \alpha_{i-1}^*$

• $\beta_{i-1} := \delta_i^2 \beta_{i-1}^*$

• $a_i := 3\xi_i - 2\alpha_{i-1} - \frac{1}{2}\beta_{i-1}$

• $b_i := 6\xi_i - 6\alpha_{i-1} - 2\beta_{i-1}$

• $\alpha_i^* := \exp(-\bar{\xi}_i) a_i$

• $\beta_i^* := \exp(-\bar{\xi}_i) b_i - (D_{-\bar{\xi}_i} \exp)(a_i) a_i$

END

2. Initialization:

• $\alpha_0^* := v_0, \beta_0^* := \dot{v}_0, \gamma_0^* := \ddot{v}_0$

• $T_0 := 1$

4. Evaluation at $t \in [t_i - t_{i-1}],$ with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:

$$ \xi_i^{[3]}(\tau_i) = \tau_i^3 \bar{\xi}_i + (\tau_i - \tau_i^3) \alpha_{i-1} + \frac{1}{2} (\tau_i^2 - \tau_i^3) \beta_{i-1} $$

$$ h_i^{[3]}(t) = h_{i-1} \exp \xi_i^{[3]}(\tau_i) $$

Algorithm 2: 4th-Order C³ continuos POE-Spline

1. Input:

• Knot points $t_0, t_1, \dots, t_N$

• Samples $h_0, h_1, \dots, h_n$

• Initial values $v_0, \dot{v}_0, \ddot{v}_0$

2. Initialization:

• $\alpha_0^* := v_0, \beta_0^* := \dot{v}_0, \gamma_0^* := \ddot{v}_0$

• $T_0 := 1$

3. Computation of Spline Coefficients:

FOR $i = 1, \dots, N$ DO

• $T_i := t_i - t_{i-1}$

• $\delta_i := T_i/T_{i-1}$

• $\bar{\xi}_i := \log(h_{i-1}^{-1}h_i)$

• $\alpha_{i-1} := \delta_i \alpha_{i-1}^*$

• $\beta_{i-1} := \delta_i^2 \beta_{i-1}^*$

• $\gamma_{i-1} := \delta_i^3 \gamma_{i-1}^*$

• $a_i := 4\xi_i - 3\alpha_{i-1} - \beta_{i-1} - \frac{1}{12} [\alpha_{i-1}, \beta_{i-1}] - \frac{1}{6}\gamma_{i-1}$

• $b_i := 12\xi_i - 12\alpha_{i-1} - 5\beta_{i-1} - \frac{1}{2} [\alpha_{i-1}, \beta_{i-1}] - \gamma_{i-1}$

• $c_i := 24\xi_i - 24\alpha_{i-1} - 12\beta_{i-1} - \frac{3}{2} [\alpha_{i-1}, \beta_{i-1}] - 3\gamma_{i-1}$

• $\alpha_i^* := \exp(-\bar{\xi}_i) a_i$

• $\beta_i^* := \exp(-\bar{\xi}_i) b_i - (D_{-\bar{\xi}_i} \exp)(a_i) a_i$

• $\gamma_i^* := \exp(-\bar{\xi}_i) c_i - 2(D_{-\bar{\xi}_i} \exp)(a_i) b_i - (D_{-\bar{\xi}_i} \exp)(b_i) a_i + (D_{-\bar{\xi}_i}^2 \exp)(a_i) (a_i) a_i$

END

4. Evaluation at $t \in [t_i - t_{i-1}],$ with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:

$$ \xi_i^{[4]}(\tau_i) = \tau_i^4 \bar{\xi}_i + (\tau_i - \tau_i^4) \alpha_{i-1} + \frac{1}{2} (\tau_i^2 - \tau_i^4) \beta_{i-1} + \frac{1}{6} (\tau_i^3 - \tau_i^4) (\gamma_{i-1} + \frac{1}{2} [\alpha_{i-1}, \beta_{i-1}]) $$

$$ h_i^{[4]}(t) = h_{i-1} \exp \xi_i^{[4]}(\tau_i) $$
## 4.2 | A $k$th-Order $C^{k-2}$ Continuos POE-Spline with given Velocities at Interpolation Points

The POE-spline in section 4.1 allows prescribing the initial velocity and derivatives. It does not admit prescribing terminal values. It is often desirable, however, to prescribe the terminal velocity, and moreover to interpolate through given states, i.e. given $h$ and $v$ at the sample points. This is addressed in the following.

In addition to $h_0, h_1, \dots, h_n \in G$, now the velocities $v_0, v_1, \dots, v_n \in g$ at the interpolation points are given. A POE-spline is derived by concatenating 2-point boundary value interpolations as they allow for prescribing $v$ and its derivatives at the start and terminal end of a segment, i.e. at the knot points. The 3rd- and 4th-order POE-splines are presented here. The general $k$th-order case becomes rather involved, and seems not to have practical relevance.

The 3rd-order interpolation (27) is used within a segment. In segment $i$ with $\tau \in [t_{i-1}, t_i]$, the initial velocity is $v_{i-1}$ and terminal velocity is $v_i$, replacing $v_0$ and $v_1$ in (27). Therewith, the $C^1$ continuity of the curve obtained by concatenating $n$ such interpolations is ensured. The corresponding computational scheme is shown as Algorithm 3. A 4th-order spline is constructed using (29) to interpolate in segment $i$. The initial and terminal values in (29) are set to $v_{i-1}$ and $v_i$ at the knot points of segment $i$. The derivative $v'_{i-1}$ at the start of segment $i$ is again calculated from the interpolation in segment $i-1$ using (35). The so obtained spline curve is $C^2$ continuos. This gives rise to Algorithm 4.

<caption>**Algorithm 3: 3rd-Order $C^1$ continuos POE Spline with given Velocities**</caption>

**1. Input:**

* Knot points $t_0, t_1, \dots, t_N$

* Samples $h_0, h_1, \dots, h_n$

* Velocities $v_0, v_1, \dots, v_n$

**3. Computation of Spline Coefficients:**
**FOR $i = 1, \dots, N$ DO**

* $T_i := t_i - t_{i-1}$

* $\delta_i := T_i/T_{i-1}$

* $\tilde{\xi}_i := \log(h_{i-1}^{-1} h_i)$

* $\alpha_i^* := T_i V_i$

* $\alpha_{i-1} := T_i V_{i-1}$

* $a_i := \exp(-\tilde{\xi}_i^{-1} \alpha_i^*)$

* $d_i := \exp(-\xi_{i-1}^{-1} \alpha_{i-1})$

**2. Initialization:**

* $\alpha_0^* := v_0$

* $T_0 := 1$

**4. Evaluation at $t \in [t_i - t_{i-1}]$, with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:**

$$ \xi_i^{[3]}(\tau_i) = (3\tau_i^2 - 2\tau_i^3) \tilde{\xi}_i + \tau_i(\tau_i - 1)^2 d_i + (\tau_i^3 - \tau_i^4) a_i \\ h_i^{[3]}(t) = h_{i-1} \exp \xi_i^{[3]}(\tau_i) $$

<caption>**Algorithm 4: 4th-Order $C^2$ continuos POE Spline with given Velocities**</caption>

**1. Input:**

* Knot points $t_0, t_1, \dots, t_N$

* Samples $h_0, h_1, \dots, h_n$

* Velocities $v_0, v_1, \dots, v_n$

* Initial derivative $\dot{v}_0$

**3. Computation of Spline Coefficients:**
**FOR $i = 1, \dots, N$ DO**

* $T_i := t_i - t_{i-1}$

* $\delta_i := T_i/T_{i-1}$

* $\tilde{\xi}_i := \log(h_{i-1}^{-1} h_i)$

* $\alpha_i^* := T_i V_i$

* $\alpha_{i-1} := T_i V_{i-1}$

* $\beta_{i-1} := \delta_i^2 V_{i-1}$

* $a_i := \exp(-\tilde{\xi}_i^{-1} \alpha_i^*)$

* $b_i := -12\tilde{\xi}_i + 6\alpha_{i-1} + \beta_{i-1} + 6a_i$

* $\beta_i^* := A_1(\tilde{\xi}_i) b_i + A_2(\tilde{\xi}_i, a_i) a_i$

END

**2. Initialization:**

* $\alpha_0^* := v_0$

* $\beta_0^* := v'_0$

* $T_0 := 1$

**4. Evaluation at $t \in [t_i - t_{i-1}]$, with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:**

$$ \xi_i^{[4]}(\tau_i) = (4\tau_i^3 - 3\tau_i^4) \tilde{\xi}_i + (1 - \tau_i)^2 (\tau_i + 2\tau_i^2) \alpha_{i-1} + \frac{1}{2}\tau_i^2 (\tau_i - 1)^2 \beta_{i-1} + (\tau_i^4 - \tau_i^5) a_i \\ h_i^{[4]}(t) = h_{i-1} \exp \xi_i^{[4]}(\tau_i) $$

**Example 2.** Interpolation of the spatial rigid body motion through a set of four given poses $h_0, \dots, h_3$ with desired velocities $v_0, \dots, v_3$, within a duration of $T = 1$ s, is considered. The rigid body poses are represented as elements of the direct product Lie group $SO(3) \times \mathbb{R}^3$. This group is compact, in contrast to $SE(3)$, thus its one-parameter subgroups define geodesics w.r.t. a bi-invariant metric on $SO(3)$ and $\mathbb{R}^3$, which are sometimes called double geodesics. This representation implies that rotations and translations are decoupled. A typical element is $h = (\mathbf{R}, \mathbf{r})$, and the corresponding left-invariant vector field is $v = (\omega, \dot{\mathbf{r}})$,
where $\omega \in so(3)$ is the body-fixed angular velocity defined via $\tilde{\omega} = \mathbf{R}^T \dot{\mathbf{R}} \in so(3)$, and $\dot{\mathbf{r}} \in \mathbb{R}^3$ is the translational velocity of the body-fixed reference frame, and $\mathbf{v}$ is called the twist in mixed representation [48]. Multiplication on this group is defined as $h_1h_2 = (\mathbf{R}_1, \mathbf{r}_1) \cdot (\mathbf{R}_2, \mathbf{r}_2) = (\mathbf{R}_1\mathbf{R}_2, \mathbf{r}_1 + \mathbf{r}_2)$. The right-trivialized differential on $SO(3)$ admits the closed form (A2). Pose $h_i$ is represented by a frame $\mathcal{F}_i$, as shown in Fig. 4. They are defined by position vectors $\mathbf{r}_i$ (represented in $\mathcal{F}_0$), and the rotation matrices $\mathbf{R}_{0i}$ (transforming from $\mathcal{F}_i$ to $\mathcal{F}_0$) as follows (all values in SI units)

$$ \mathbf{r}_1 = \begin{bmatrix} 1 \\ 4 \end{bmatrix}, \mathbf{r}_2 = \begin{bmatrix} 4 \\ 4 \end{bmatrix}, \mathbf{r}_3 = \begin{bmatrix} 8 \\ 4 \end{bmatrix}, \mathbf{R}_{01} = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \mathbf{R}_{02} = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}, \mathbf{R}_{03} = \mathbf{I}. $$

The angular velocity is required to be zero at all four points, $\omega_i = 0, i = 0, \dots, 3$. The initial and terminal linear velocity is zero, $\dot{\mathbf{r}}_0 = \dot{\mathbf{r}}_3 = \mathbf{0}$, and at the via points is prescribed as $\dot{\mathbf{r}}_1 = [10, 0, 0]^T$ and $\dot{\mathbf{r}}_2 = [0, 0, 10]^T$, and the initial velocity as $\dot{\mathbf{r}}_0 = \mathbf{0}$. The results obtained with the 3rd- and 4th-order spline are shown in Fig. 4. For the 4th-order spline, the derivative $v_0'$ is prescribed with zero angular acceleration $\dot{\omega}_0 = \mathbf{0}$ and linear acceleration $\ddot{r}_0 = [-10, 0, 0]^T$. Shown is the moving frame at equidistant samples of $t \in [0, T]$. The difference of the 3rd- and 4th-order curve are due to fact that the 4th-order curve respects continuity of the jerk $v''$. The special situation when the velocity at all points is set to zero is shown in Fig. 5. Then, the 3rd-order spline interpolation reduces to the concatenation of 3rd-order 2-point interpolations with matching velocities. When using the 4th-order interpolations, additionally the accelerations match at the via points, explaining the apparent difference of the curves. For comparison, the motion curve constructed with a De Casteljau algorithm on $SO(3) \times \mathbb{R}^3$ is shown in Fig. 6. The attitude, i.e. curve in $SO(3)$, is constructed with algorithm 9 (appendix B), and the position, i.e. curve in vector space $\mathbb{R}^3$, is constructed with the classical De Casteljau algorithm. It would require a Bézier-spline algorithm to construct a motion passing through the two intermediate configurations, which is still topic of research [30, 57].

<figure><img src="image_5.png" /><figcaption>FIGURE 4 Interpolation of a spatial rigid body motion through four given frame locations $\mathcal{F}_0, \dots, \mathcal{F}_3$, with prescribed velocities. a) Results when using 3rd-order, and b) when using the 4th-order spline. $\mathcal{F}_0$ is the global reference frame.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 5 Interpolation of a spatial rigid body motion through four given frame locations $\mathcal{F}_0, \dots, \mathcal{F}_3$ when the velocity is zero at the given frame locations, using a) 3rd-order, and b) 4th-order spline. $\mathcal{F}_0$ is the global reference frame.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6 Bézier curve interpolating between start and end pose, $\mathcal{F}_0$ and $\mathcal{F}_3$, of the spatial rigid body motion constructed with the De Casteljau algorithm on $SO(3) \times \mathbb{R}^3$ using intermediate frames $\mathcal{F}_1$ and $\mathcal{F}_2$ as control points.</figcaption></figure>

## 5 | A LIMITATION OF POE-SPLINES

A limitation inherent to $k$th-order POE-spline with zero initial local coordinates $\xi_i^{[k]}(0) = 0$ is that they cannot exactly reconstruct a $k$th-order curve on $G$, i.e. a curve $h(t) = h_0 \exp \xi(t)$ where $\xi(t)$ is a $k$th-order polynomial in $t$. Assume that the curve is to be interpolated through $h_i = h(t_i), i = 0, \dots, n$. In interval $t \in [t_{i-1}, t_i]$, the curve can be expressed as $h(t) = h_0 \exp \xi(t) = h_0 \exp \xi(t_{i-1}) \exp \eta(t) = h_{i-1} \exp \eta(t)$. The coordinates $\xi(t)$ and $\eta(t)$ are related by the BCH formula as

$$ \eta(t) = \xi(t) - \xi(t_{i-1}) - \frac{1}{2} [\xi(t_{i-1}), \xi(t)] + \frac{1}{12} [\xi(t_{i-1}), [\xi(t_{i-1}), \xi(t)]] - \frac{1}{12} [\xi(t), [\xi(t), \xi(t_{i-1})]] + \dots \quad (40) $$

If $\xi(t)$ is a $k$th-order polynomial in $t$, then $\eta(t)$ is generally not a $k$th-order polynomial, which requires all nested Lie brackets in (40) to vanish. In case of standard spline interpolation on vector spaces, and on general Abelian groups, this condition is satisfied, of course.

**Proposition 2.** A $k$th-order POE-spline can only exactly reconstruct a curve of degree $k$ defined on a $1$-parameter subgroup of $G$.

This limitation has not been addressed in the literature. It is particularly relevant when a curve in $G$ is to be interpolated, rather than to be generated from given via point. Besides this issue, another limitation is that the initial pose $h_0$ must be known beforehand, rather than being the outcome of the interpolation.
**Example 3.** To demonstrate proposition 2, and thus lemma 2, the interpolation of a spatial rotation in SO (3) is considered. The curve to be interpolated in SO(3) is expressed as $h(t) = \exp(\xi(t))$, where $h(t) \in SO(3)$ is the rotation matrix, and $\xi(t)$ the instantaneous rotation axis times angle, and the exponential map is the Euler-Rodrigues formula (A1). First consider a motion in the 1-parameter subgroup of rations about the constant (for simplicity non-unit) axis $\xi_0 = [0.5, 1.5, 1] \in \mathbb{R}^3$. The particular motion to be interpolated is described by $\xi(t) = (t^3 - t^2 + 3t)\xi_0$, which is cubic in $t$. As input to the spline interpolation, $n = 11$ equidistant samples $h_0, \dots, h_{10} \in SO(3)$ are computed along the curve at $\tau_i, i = 0, \dots, 1$ giving rise to a 10-segment spline. The 3rd-order cubic POE-spline with given initial values (Algorithm 1) and with given velocities at the sample points (Algorithm 3) are applied. For the latter, also the corresponding angular velocities $v_0, \dots, v_{10} \in SO(3)$ are computed as $v_i = h^{-1}(\tau_i) \dot{h}(\tau_i)$. The interpolation error is computed as $\epsilon = \|\log(\Delta h)\|$, with $\Delta h = h^{-1}(t) h^{[3]}(t)$, where $h^{[3]}(t) = \exp(\xi^{[3]}(t))$ is the curve obtained with the 3rd-order interpolation. Fig. 7a) shows the result obtained with the two 3rd-order interpolations. Clearly, this 3rd-order motion can be reconstructed (numerically) exactly by the 3rd-order spline. Now consider the spatial rotation described by $\xi(t) = t\xi_1 + t^3\xi_2$, with $\xi_1 = [0.1, 0, 0.2] \in \mathbb{R}^3$ and $\xi_2 = [0, 1.5, 0] \in \mathbb{R}^3$. It is obvious from Fig. 7b) that the 3rd-order spline cannot reconstruct this cubic motion. Also the De Casteljau algorithm 9 (appendix B) is applied for comparison. For both motions the error magnitude is about 0.06 rad. It is therefore not shown in the figures.

<figure><img src="image_8.png" /><figcaption>FIGURE 7 Error when interpolating a) a cubic 1-parameter motion in SO(3) and b) a general cubic motion in SO(3) with a 3rd-order POE spline.</figcaption></figure>

## 6 | 2-POINT INTERPOLATION WITH NON-ZERO INITIAL VALUE $\xi(0) \neq 0$

If $\xi(0)$ is non-zero, the series solution (7) does not apply. The following generalized the series solution, with (7) as special case.

**Theorem 1.** i) The solution of the left-invariant normalized Poisson equation (3) can locally be expressed as $g(\tau) = g_0 \exp \xi(\tau)$, where $\xi(\tau)$ satisfies the ODE (5) with $\xi(0) \neq 0$ small. ii) At $\tau = 0$, the solution admits the series expansion

$$ \xi(\tau) = \sum_{k \ge 0} \frac{\tau^k}{k!} \Phi_k(\xi(0), v(0), v'(0), \dots, v^{(k-1)}(0)) \quad (41) $$

with $\Phi_0 = \xi_0$, and $\Phi_k(\xi(\tau), v(\tau), v'(\tau), \dots, v^{(s-1)}(\tau))$ in (23). iii) If the initial value is zero, $\xi(0) = 0$, the coefficients are given in (8), i.e. $\Phi_k = \Phi_k^0, k \ge 1$ and $\Phi_k = 0$.

*Proof.* i) With $g(\tau) = g_0 \exp \xi(\tau)$, the left Poisson equation (3) becomes $v(\tau) = \exp(-\xi(\tau)) \frac{d}{d\tau} \exp(-\xi(\tau)) = \mathrm{dexp}_{-\xi(\tau)} \xi'(\tau)$. This relation is independent of whether $\xi(0) = 0$. The dexp map exists for $\xi \in g$ small. ii) Inserting the series expansion of $\xi'(\tau)$ and $\mathrm{dexp}_{-\xi(\tau)} v(\tau)$ at $\tau = 0$ in (5) yields $\sum_{j \ge 0} \frac{\tau^j}{j!} \Phi_{j+1} = \mathrm{dexp}_{-\xi(0)} v(0) + \sum_{j \ge 1} \frac{\tau^j}{j!} \frac{d^j}{d\tau^j} \left( \mathrm{dexp}_{-\xi(\tau)} v(\tau) \right) \Big|_{\tau=0}$. Application of Leibnitz rule to the last term, and comparing coefficients for powers of $\tau$ yields (23) and (24). The dexp map being analytic in a neighborhood of the identity ensures convergence. Convergence proofs can be found in [44, 70]. iii) This was proven in [47]. □
**Remark 7.** The condition that $\xi(0) \neq 0$ is small ensures that $\exp \xi(\tau)$ remains in a neighborhood of the identity in $G$. In practical applications this requirement is lifted, and parameter singularities are tackled using the limit values of the respective exp map.

**Remark 8.** It is instructive to consider case iii) as the limit of $\Phi_k$ in (23) for $\xi \to 0$. To this end, the series expansion (6) is used. The series expansions of the first and second derivatives are, for instance,

$$ (D_{\xi} \exp^{-1})(\eta) = \sum_{i \ge 1} \frac{B_i}{i!} \sum_{\substack{j,l \ge 0 \\ j+l+1=i}} ad_{\xi}^{j} ad_{\eta}^{l} ad_{\xi}^{l} \quad (42) $$

$$ (D_{\xi}^2 \exp^{-1})(\eta)(\omega) = \sum_{i \ge 2} \frac{B_i}{i!} \sum_{\substack{j,l \ge 0 \\ j+l+2=i}} \left( \sum_{\substack{r,s \ge 0 \\ r+s=j}} ad_{\xi}^r ad_{\omega}^s ad_{\xi}^l ad_{\eta}^l ad_{\xi}^r ad_{\eta}^r ad_{\omega}^s ad_{\xi}^s + \sum_{\substack{r,s \ge 0 \\ r+s=l}} ad_{\xi}^j ad_{\eta}^l ad_{\xi}^r ad_{\eta}^l ad_{\omega}^s ad_{\xi}^s \right) \quad (43) $$

with $\xi, \eta, \omega \in g$. The third derivative vanishes identically since $B_3 = 0$. Therewith, the limits at $\xi(0) = 0$ are found as

$$ (D_0 \exp^{-1})(\eta) = -\frac{1}{2} ad_{\eta}, \quad (D_0^2 \exp^{-1})(\eta)(\omega) = \frac{1}{12} (ad_{\omega} ad_{\eta} + ad_{\eta} ad_{\omega}). \quad (44) $$

The coefficients (10) are then obtained from (44).

Theorem 1 gives rise to kth-order 2-point interpolation formulae with $\xi(0) \neq 0$. The derivation is analogous to that in section 3. Denote $\Delta\xi := \xi - \xi_0$, the initial value interpolation (16) becomes

$$ \xi^{[k]}(\tau) = \xi_0 + \tau^k \Delta\xi - \sum_{j=1}^{k-1} \frac{1}{j!} (\tau^j - \tau^k) \mathbf{v}_0^{(j-1)} + \sum_{r=3}^{k} (\tau^r - \tau^k) \mathbf{a}_r \quad (45) $$

and in particular, the 3rd- and 4th-order interpolation formulae are

$$ \xi^{[3]}(\tau) = \xi_0 + \tau^3 \Delta\xi + (\tau - \tau^3) \mathbf{v}_0 + \frac{1}{2} (\tau^2 - \tau^3) \mathbf{v}_0' \quad (46) $$

$$ \xi^{[4]}(\tau) = \xi_0 + \tau^4 \Delta\xi + (\tau - \tau^4) \mathbf{v}_0 + \frac{1}{2}(\tau^2 - \tau^4) \mathbf{v}_0' + \frac{1}{6}(\tau^3 - \tau^4) (\mathbf{v}_0'' + \frac{1}{2}[\mathbf{v}_0', \mathbf{v}_0]). $$

The 3rd- and 4th-order boundary value interpolation, (27) and (29) respectively, become

$$ \xi^{[3]}(\tau) = \xi_0 + (3\tau^2 - 2\tau^3) \Delta\xi + \tau(\tau - 1)^2 \mathbf{v}_0 + (\tau^3 - \tau^2) \exp(-\xi) \mathbf{v}_1 \quad (47) $$

$$ \xi^{[4]}(\tau) = \xi_0 + (4\tau^3 - 3\tau^4) \Delta\xi + \tau(\tau - \tau^2)(1 + 2\tau) \mathbf{v}_0 + (\tau^4 - \tau^3) \exp(-\xi) \mathbf{v}_1 + \frac{1}{2}\tau^2(1 - \tau)^2 \mathbf{v}_0'. \quad (48) $$

# 7 | GLOBALLY PARAMETERIZED KTH-ORDER CONTINUOS SPLINES

Instead of concatenating local interpolations by group multiplication, which leads to the POE-spline, the interpolation is defined by a single exponential of a curve in $g$, represented by $\xi(t)$ defined on the entire range $t \in [0, T]$ of the global path parameter.

**Definition 3.** A kth-order $C^p$-continuos global spline interpolating through $h_0, h_1, \dots, h_n \in G$ is a map

$$ \bar{h}(t) = \bar{h}_0 \exp \xi(t), \quad t \in [0, T] \quad (49) $$

satisfying the matching condition $h(t_i) = h_i$ and the initial conditions $\bar{v}^{(s)}(0) = \mathbf{v}_0^{(s)}$, $s = 0, \dots, p-1$, where the coordinate function $\xi(t)$ and $\bar{v}^{(s)}(t)$, $s = 0, \dots, p-1$ are continuos, with left-invariant vector field defined as $\bar{v}(t) := \bar{h}(t)^{-1} \bar{h}'(t)$.

This yields a global interpolation in terms of local coordinates $\xi(t)$ on $G$. The canonical coordinates will still be interpolated separately for each spline segment, but the propagation through the segment is described as $\xi(t) = \xi(t_{i-1}) + \xi_i^{[k]}(\tau_i)$, $t \in [t_{i-1}, t_i]$, which does not suffer from the issue discussed in section 5 as the concatenation is done on the Lie algebra (a vector space). Notice that global refers to fact that local coordinates $\xi$ are used to parameterize the complete curve in $G$. Such global interpolations
were addressed e.g. in [27] for rigid body motions. Since $\xi(0)$ is no longer required to be zero, in general $\bar{h}_0 \neq h_0$. Nevertheless, $\xi(0) = 0$ is an obvious choice.

## 7.1 | A $k$th-Order $C^{k-1}$ Continuos Spline with given Initial Values

The key to this spline formulation is the 2-point interpolation in section 6. The POE-formulation is readily amended leading to Algorithms 5 and 6. The canonical coordinates at the end of segment $i$ are now obtained from $h_0$ and $h_i$ as $\tilde{\xi}_i := \log(h_0^{-1}h_i)$, rather than from the relative local increment $h_0^{-1}h_i$ as in case of POE splines.

*Remark 9.* The uniqueness of the log map on the Lie group $G$ must be taken into account in numerical computations when the 'distance' of the interpolation points $h_i$ becomes large. In case of $SO(3)$, and thus $SE(3)$, the log map only has a (removable) singularity at the origin, but it is only unique for $\|x\| \le \pi$. This non-uniqueness issue when extracting rotation axis and angle from a given rotation matrix (or unit quaternions) is well-known [3, 10, 72], and can be treated if continuity of the curve $h(t)$ is ensured.

**Algorithm 5: 3rd-Order $C^2$ Continuos Global Spline**

**1. Input:**

* Knot points $t_0, t_1, \dots, t_N$

* Samples $h_0, h_1, \dots, h_n$

* Initial values $\mathbf{v}_0, \mathbf{v}'_0$

* Reference $\bar{h}_0$

**3. Computation of Spline Coefficients:**

FOR $i = 1, \dots, N$ DO

* $T_i := t_i - t_{i-1}$

* $\delta_i := T_i/T_{i-1}$

* $\xi_i := \log(h_0^{-1}h_i)$

* $\Delta\xi_i := \xi_i - \xi_{i-1}$

* $\alpha_{i-1} := \delta_i\alpha_{i-1}^*$

* $\beta_{i-1} := \delta_i^2\beta_{i-1}^*$

* $\mathbf{d}_i := \exp_{-\xi_{i-1}} \alpha_{i-1}$

* $\mathbf{e}_i := \exp_{-\xi_{i-1}} \beta_{i-1} - (D_{-\xi_{i-1}} \exp^{-1})(\mathbf{d}_i) \alpha_{i-1}$

* $\mathbf{a}_i := 3\xi_i - 2\alpha_{i-1} - \frac{1}{2}\beta_{i-1}$

* $\mathbf{b}_i := 6\xi_i - 6\alpha_{i-1} - 2\beta_{i-1}$

* $\alpha_i^* := \exp_{-\xi_i} \mathbf{a}_i$

* $\beta_i^* := \exp_{-\xi_i} \mathbf{b}_i - (D_{-\xi_i} \exp^{-1})(\mathbf{a}_i) \mathbf{a}_i$

END

**2. Initialization:**

* $\alpha_0^* := \mathbf{v}_0, \beta_0^* := \mathbf{v}'_0, \gamma_0^* := \mathbf{v}''_0$

* $\beta_0^* := \mathbf{v}'_0$

* $T_0 := 1$

**4. Evaluation at $t \in [t_{i-1}, t_i]$, with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:**

$$ \xi_i^{[3]}(\tau_i) = \xi_{i-1} + \tau_i^3 \Delta\xi_i + (\tau_i - \tau_i^3) \mathbf{d}_i + \frac{1}{2} (\tau_i^2 - \tau_i^3) \mathbf{e}_i \\ h^{[3]}(t) = \bar{h}_0 \exp \xi_i^{[3]}(\tau_i(t)) $$
**Algorithm 6: 4th-Order C³ Continuos Global Spline**

**1. Input:**

* Knot points $t_0, t_1, \dots, t_N$

* Samples $h_0, h_1, \dots, h_n$

* Initial values $v_0, v'_0, v''_0$

* Reference $\bar{h}_0$

**2. Initialization:**

* $\alpha_0^* := v_0, \beta_0^* := v'_0, \gamma_0^* := v''_0$

* $\beta_0^* := v'_0$

* $\gamma_0^* := v''_0$

* $T_0 := 1$

**3. Computation of Spline Coefficients:**

FOR $i = 1, \dots, N$ DO

* $T_i := t_i - t_{i-1}$

* $\delta_i := T_i/T_{i-1}$

* $\xi_i := \log(h_0^{-1}h_i)$

* $\Delta\xi_i := \xi_i - \xi_{i-1}$

* $\alpha_{i-1} := \delta_i \alpha_{i-1}^*$

* $\beta_{i-1} := \delta_i^2 \beta_{i-1}^*$

* $\gamma_{i-1} := \delta_i^3 \gamma_{i-1}^*$

* $d_i := \exp(-\xi_{i-1}) \alpha_{i-1}$

* $e_i := \exp(-\xi_{i-1}) \beta_{i-1} - (D_{-\xi_{i-1}} \exp^{-1})(d_i) \alpha_{i-1}$

* $f_i := \exp(-\xi_{i-1}) \gamma_{i-1} - 2(D_{-\xi_{i-1}} \exp^{-1})(d_i) \beta_{i-1} - (D_{-\xi_{i-1}} \exp^{-1})(e_i) \alpha_{i-1} + (D_{-\xi_{i-1}}^2 \exp)(d_i, d_i) \alpha_{i-1}$

* $a_i := 4\Delta\xi_i - 3d_i - e_i - \frac{1}{6}f_i$

* $b_i := 12\Delta\xi_i - 12d_i - 5e_i - f_i$

* $c_i := 24\Delta\xi_i - 24d_i - 12e_i - 3f_i$

* $\alpha_i^* := \exp(-\xi_i) a_i$

* $\beta_i^* := \exp(-\xi_i) b_i - (D_{-\xi_i} \exp)(a_i) a_i$

* $\gamma_i^* := \exp(-\xi_i) c_i - 2(D_{-\xi_i} \exp)(a_i) b_i - (D_{-\xi_i} \exp)(b_i) a_i + (D_{-\xi_i}^2 \exp)(a_i, a_i) a_i$

END

**4. Evaluation** at $t \in [t_{i-1}, t_i]$, with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:

$$ \xi_i^{[4]}(\tau_i) = \xi_{i-1} + \tau_i^4 \Delta\xi_i + (\tau_i - \tau_i^4) d_i + \frac{1}{2}(\tau_i^2 - \tau_i^4) e_i + \frac{1}{6}(\tau_i^3 - \tau_i^4) f_i $$

$$ \bar{h}_i^{[4]}(t) = \bar{h}_0 \exp \xi_i^{[4]}(\tau_i(t)) $$

## 7.2 | A kth-Order C$^{k-2}$ Continuos Spline with given Velocities at Interpolation Points

As in section 4.2, it is now assumed that velocities $v_0, v_1, \dots, v_n$ are prescribed at the knot points. Algorithms 3 and 4 are readily adopted to the case of non-zero initial value $\xi_i(0) \neq 0$. This leads to Algorithms 7 and 8. The splines are constructed by concatenating the 2-point boundary interpolation formulae (27) and (29). In the 4th-order algorithm, the derivative $v'_i(t_{i-1})$ at begin of segment i is computed with (38).
Algorithm 7: 3rd-Order C¹ Continuos Global Spline with Given Velocities

1. Input:

* Knot points $t_0, t_1, \dots, t_N$

* Samples $h_0, h_1, \dots, h_n$

* Velocities $\mathbf{v}_0, \mathbf{v}_1, \dots, \mathbf{v}_n$

* Reference $\bar{h}_0$

2. Initialization:

* $\alpha_0^* := \mathbf{v}_0$

* $T_0 := 1$

3. Computation of Spline Coefficients:

FOR $i = 1, \dots, N$ DO

* $T_i := t_i - t_{i-1}$

* $\delta_i := T_i/T_{i-1}$

* $\xi_i := \log(h_0^{-1}h_i)$

* $\Delta\xi_i = \log(h_{i-1}^{-1}h_i)$

* $\alpha_i^* := T_i V_i$

* $\alpha_{i-1} := T_i V_{i-1}$

* $\alpha_i := \exp(-\xi_i^{-1} \alpha_i^*)$

* $\mathbf{d}_i := \exp(-\xi_{i-1}^{-1} \alpha_{i-1})$

END

4. Evaluation at $t \in [t_{i-1}, t_i]$, with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:

$$ \xi_i^{[3]}(\tau_i) = \bar{\xi}_i + \tau_i^3 \Delta\xi_i + (\tau_i - \tau_i^3) \alpha_{i-1} + \frac{1}{2} (\tau_i^2 - \tau_i^3) \beta_{i-1} $$

$$ h_i^{[3]}(t) = \bar{h}_0 \exp \xi_i^{[3]}(\tau_i(t)) $$

Algorithm 8: 4th-Order C² Continuos Global Spline with given Velocities

1. Input:

* Knot points $t_0, t_1, \dots, t_N$

* Samples $h_0, h_1, \dots, h_n$

* Velocities $\mathbf{v}_0, \mathbf{v}_1, \dots, \mathbf{v}_n$

* Initial derivative $\mathbf{v}_0'$

* Reference $\bar{h}_0$

2. Initialization:

* $\alpha_0^* := \mathbf{v}_0$

* $\beta_0^* := \mathbf{v}_0'$

* $T_0 := 1$

3. Computation of Spline Coefficients:

FOR $i = 1, \dots, N$ DO

* $T_i := t_i - t_{i-1}$

* $\delta_i := T_i/T_{i-1}$

* $\bar{\xi}_i := \log(h_0^{-1}h_i)$

* $\Delta\bar{\xi}_i = \log(h_{i-1}^{-1}h_i)$

* $\alpha_i^* := T_i V_i$

* $\alpha_{i-1} := T_i V_{i-1}$

* $\alpha_i := \exp(-\xi_i^{-1} \alpha_i^*)$

* $\mathbf{b}_i := -12\xi_i + 6\alpha_{i-1} + \beta_{i-1} + 6\mathbf{a}_i$

* $\beta_i^* := \exp(-\xi_i \mathbf{b}_i - (\mathbf{D}_i \exp \xi_i)(\mathbf{a}_i) \mathbf{a}_i)$

END

4. Evaluation at $t \in [t_{i-1}, t_i]$, with $\tau_i(t) = \frac{t-t_{i-1}}{T_i}$:

$$ \xi_i^{[4]}(\tau_i) = \xi_{i-1} + (4\tau_i^3 - 3\tau_i^4) \Delta\xi_i + (1 - \tau_i)^2 (\tau_i + 2\tau_i^2) \mathbf{d}_i + \frac{1}{2}\tau_i^2 (\tau_i - 1)^2 \mathbf{e}_i + (\tau_i^4 - \tau_i^3) \mathbf{a}_i $$

$$ h_i^{[4]}(t) = \bar{h}_0 \exp \xi_i^{[4]}(\tau_i(t)) $$

**Example 4.** The interpolation of a cubic spatial rotation in SO (3), described by $\xi(t) = t\xi_1 + t^3\xi_2$, was investigated in example 3, revealing the limitation of POE-splines. The same input data, $h_0, \dots, h_{10}$ and $\mathbf{v}_0, \dots, \mathbf{v}_{10}$, are now used for a global 3rd-order spline interpolation. The results obtained with the 3rd-order global spline for given initial values (Algorithm 5) and when velocities samples are provided (Algorithm 7) are shown in Fig. 8. The global interpolation can clearly reconstruct the cubic curve, in contrast to the POE-spline (see Fig. 7). The same is true when applying the 4th-order interpolation.

**Example 5.** The shape reconstruction of a slender rod undergoing spatial deflections is addressed. In example 1, the shape was approximated by a 2-point interpolation method. Now a spline on SE (3) is used to represent the rod, and parameterized by the normalized global path parameter $t \in [0, 1]$. To this end, the cross section displacement is sampled at $n = 5$ equidistant points $t_1, \dots, t_n$ along the rod yielding $h_0, \dots, h_n \in SE(3)$. The poses are obtained from the reference solution $\xi(t)$ as $h_i = \exp \xi(t_i)$ (see example 1). For this example, the deformation to be interpolated is the result of a concentrated load $\mathbf{W} = [M, F]$ in (32), applied at the terminal end of the rod ($\tau = 1$). Deformations due to a pure moment about a principal cross section axis are exactly recovered by any of the splines since the deformation is described by a constant screw axis (constant curvature). In the following,
<figure><img src="image_9.png" /><figcaption>FIGURE 8 Error when interpolating a general cubic motion in SO(3) with a global 3rd-order spline.</figcaption></figure>

the applied wrench is computed with the force **F** (1) = [0, 0.12, -0.08] and moment **M** (1) = [0, -0.008, -0.012]. The resulting shape is shown in Fig. 9.

The initial and terminal deformation measure, $v_0$ and $v_1$, and the derivative $v'_0$, needed for the spline interpolation, are computed from the reference solution. Fig. 10a) shows the interpolation errors when using the 3rd- and 4th-order POE spline algorithms 1 and 2, respectively. The global spline, constructed with algorithms 5 and 6, yields the results in Fig. 10b). Clearly, for this general deformation the POE- and the global spline yield almost the same results. The visible larger error of the 4th-order spline can be explained by its sensitivity to the derivative $v'_0$, which is obtained numerically. For comparison, the result of a Bézier curve interpolation constructed with the De Casteljau algorithm 9 (appendix B) is shown in blue. The 5 samples yield a 4th-order Bézier curve. Apparently, the curve does not pass through the intermediate configurations, which serve as control points, but interpolates between start and end configuration (see also appendix B). To make the curve pass through the points needs computation of control points in the Lie group, which is an open problem and not the topic of this paper. For constructing B-splines on SE (3), a De Boor algorithm was presented in [21, 67]. Next the POE spline algorithms 3 and 4 are applied. To this end, the deformation $v_i$, $i = 0, ..., n$ is computed at the $n = 5$ interpolation points. This is found from the equilibrium conditions as $v_i = K^{-1}Ad_{h_i}^T K(v_0 - v^0) + v^0$, where $Ad_{h_i}$ is the adjoined operator on SE (3) [53, 62]. The improvement in accuracy when providing the deformations $v_i$ at the interpolation points is apparent from Fig. 11a). This alleviates the sensitivity of the 4th-order spline, which better approximates the shape. The results obtained with the global spline algorithms 5 and 6, with given deformations $v_i$, $i = 0, ..., n$, are shown in Fig. 11b), which are again similar to those of the POE-spline. Finally, Fig. 12 shows the interpolation error when $n = 10$ interpolation points are used, i.e. the number of segments is doubled.

<figure><img src="image_10.png" /><figcaption>FIGURE 9 Static deformation of the slender rod when subjected to the wrench **W** (1) = [**M** (1), **F** (1)] at the terminal end.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**FIGURE 10** Interpolation error $\epsilon_r$ and $\epsilon_p$ for orientation and position when using a) POE-spline, and b) global spline. Additionally shown in blue is the error for a Bézier curve interpolating through the 5 points using De Casteljau algorithm 9.</figcaption></figure>

While for this example the $h_i$ and $v_i$ are computed from a reference solution, it should be remarked that it depends on the application how the $h_i$ are obtained.

<figure><img src="image_12.png" /><figcaption>**FIGURE 11** Interpolation error $\epsilon_r$ and $\epsilon_p$ for orientation and position obtained with a) POE-spline, and b) global spline, when deformations are prescribed at the interpolation points.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**FIGURE 12** Interpolation error $\epsilon_r$ and $\epsilon_p$ for orientation and position obtained with a) POE-spline, and b) global spline with 20 segments, when deformations are prescribed at the interpolation points.</figcaption></figure>
**Example 6.** The elastic rod from example 5 is considered, but with a rectangular cross section 4 × 12 mm². The rod is subjected to a force **F** (0) = [0, −0.03, 0.01] and moment **M** (0) = [0, 0.002, 0.002] applied at the initial end at τ = 0 that is in static equilibrium with the force **F** (1) = [−0.0021, 0.0192, −0.0181] and moment **M** (1) = [0.0296, 0.0082, 0.0075] at the terminal with τ = 1. The numerically exact displacement obtained as solution of (3) and (32) is shown in Fig. 13 in gray color. The shape reconstructed with the 3rd-order interpolation is shown in Fig. 13a), and in Fig. 13b) when additionally the deformation is prescribed at the interpolation points. The interpolation error is shown in Fig. 14. Results are shown in Fig. 15 when the rod exhibits larger variations in the deformation field. In all examples it is apparent that providing the deformations at the interpolation points reduces the error amplification with increasing via points, which is well-known for spline interpolations in general. For completeness, the results are shown in Fig. 14 for the 4th-order Bézier curve constructed with the 5 equidistant samples along the rod used as control points using the De Casteljau algorithm 9.

<figure><img src="image_14.png" /><figcaption>**FIGURE 13** Shape of flexible rod resulting from a) the 3rd-order interpolation and b) the 3rd-order interpolation with deformation prescribed at the interpolation points (red). The exact numerical solution is shown in gray color (partially covered).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**FIGURE 14** a) Position error $\epsilon_p$ and b) orientation error $\epsilon_q$ obtained with 3r-order piecewise POE-spline (i), global POE spline (ii), piecewise POE-spline (iii) and global POE-spline (iv) with prescribed deformation at interpolation points, and the Bézier curve of order 4 constructed with the De Casteljau algorithm 9 (v).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**FIGURE 15** Shape of flexible rod with larger variations of the deformation. a) Numerically exact solution, b) 3rd-order piecewise interpolation, c) 3rd-order piecewise interpolation with prescribed deformation at the five interpolation points.</figcaption></figure>

## 8 | CONCLUSION AND OUTLOOK

This paper addressed the interpolation problem on connected Lie groups, with particular focus on interpolating rigid body motions, i.e. curves in SO(3) and SE(3). A kth-order 2-point interpolation was derived from a truncated series solution of the Poisson equation. The minimum acceleration and minimum energy curves w.r.t. a left invariant metric are recovered as spacial cases. Thereupon, kth-order $C^{k-1}$ continuos POE spline interpolation schemes can be constructed. A 3rd- and 4th-order interpolation scheme was presented. The 3rd-order spline algorithm is equivalent to that proposed in [35, 56]. An important contribution of this paper is the extension of these algorithms to additionally allow prescribing a vector field, in addition to the points, the interpolation curve must match. In case of rigid body motions, this admits prescribing the velocity at the interpolation points, thus allows trajectory planning through a given set of states rather than just poses. These formulations also reduce the sensitivity of all such spline formulations w.r.t. the initial conditions. Interpolation of frame curves in Euclidean space has various further applications. One such applications is the shape reconstruction of geometric non-linear flexible rods modeled using Cosserat theory. This leads to an interpolation problem on SE(3). The main benefit of the proposed formulation is that it does not explicitly rely on geodesics. The reported results confirm the potential for such applications.

It was pointed out that a kth-order POE spline cannot exactly represent a kth-order motion on G, i.e. a curve described by canonical coordinates that are a kth-order polynomials in the path parameter. To overcome this problem, globally parameterized splines were introduced. This rests on a 2-point interpolation with non-zero initial values of the coordinates, in contrast to existing interpolation schemes. The global parameterization is particularly relevant when a polynomial curve in G is to be interpolated, which is a rather special situation, however. The global splines are also important when the initial value of the exponential coordinates is non-zero. This happens in case of continuum elements in soft robots, for instance, where the initial $h(0) \neq h_0$ and terminal pose $h(1)$ are computed from an elastic equilibrium.

Future work will investigate the interpolation on subspaces of G that are not subgroups. An important special case are symmetric spaces [13, 73]. To this end, the higher-order solutions must be derived taking into account particular generators of the symmetric space.

## ACKNOWLEDGMENTS

This work was supported by the LCM-K2 Center within the framework of the Austrian COMET-K2 program.

## CONFLICT OF INTEREST

The authors declare no potential conflict of interests.
# APPENDIX

## A CLOSED FORM EXPRESSIONS FOR THE LIE GROUPS SO (3) AND SE (3)

In the following, the closed form relations are summarized, which can be found in [49], and partially in [53, 62]. To be consistent with the notation of this paper, Lie group elements are denoted with $g$ and elements of the Lie algebra with $\xi$.

The Lie group SO (3) is the group of orthogonal 3 × 3 matrices with determinate equal to 1, referred to as rotation matrices. A rotation matrix $g \in SO(3)$ is expressed in terms of a unit vector $n \in \mathbb{R}^3$ along the rotation axis and rotation angle $\varphi$ via the exponential map

$$ \exp \tilde{\xi} = \mathbf{I} + \alpha \tilde{\xi} + \frac{1}{2} \beta \tilde{\xi}^2 \quad (A1) $$

where $\xi := \varphi n$, $\alpha := \text{sinc}\varphi$, $\beta := \text{sinc}^2 \frac{\varphi}{2}$, and $\tilde{\xi} \in so(3)$ is the skew symmetric matrix associated with $\xi \in \mathbb{R}^3$. The Lie algebra so (3) is the vector space of skew symmetric 3 × 3 matrices. With the unique correspondence of $\xi \in \mathbb{R}^3$ and $\tilde{\xi} \in so(3)$, the notation $\xi$ and $\tilde{\xi}$ are used interchangeably. The right-trivialized differential of (A1), and its inverse admit the closed form expressions

$$ \mathrm{dexp}_{\xi} = \mathbf{I} + \frac{\beta}{2}\tilde{\xi} + \frac{1}{\|\mathbf{x}\|^2}(1-\alpha)\tilde{\xi}^2 = \mathbf{I} + \frac{\beta}{2}\tilde{\xi} + (1-\alpha)\tilde{\mathbf{n}}^2 \quad (A2) $$

$$ \mathrm{dexp}_{\xi}^{-1} = \mathbf{I} - \frac{1}{2}\tilde{\xi} + \frac{1}{\|\mathbf{x}\|^2}(1-\gamma)\tilde{\xi}^2 = \mathbf{I} - \frac{1}{2}\tilde{\xi} + (1-\gamma)\tilde{\mathbf{n}}^2 \quad (A3) $$

where $\gamma := \frac{\alpha}{\beta}$, $\delta := \frac{1-\alpha}{\varphi^2}$. The Lie bracket on so (3) is given by the cross product, $[\xi_1, \xi_2] = ad_{\xi_1}\xi_2 = \tilde{\xi}_1\xi_2 = \xi_1 \times \xi_2$, for $\xi_1, \xi_2 \in \mathbb{R}^3 \cong so(3)$.

The Lie group SE (3) is the group of isometric and orientation preserving transformations in $E^3$ [62]. In matrix representation, a typical element is

$$ g = \begin{bmatrix} \mathbf{R} & \mathbf{r} \\ \mathbf{0} & 1 \end{bmatrix} \in SE(3) \quad (A4) $$

where $\mathbf{R} \in SO(3)$, $\mathbf{r} \in \mathbb{R}^3$. The Lie algebra se (3) is the vector space of matrices

$$ \hat{\xi} = \begin{bmatrix} \tilde{\mathbf{x}} & \mathbf{y} \\ \mathbf{0} & 0 \end{bmatrix} \in se(3) \quad (A5) $$

where $\tilde{\mathbf{x}} \in so(3)$, $\mathbf{y} \in \mathbb{R}^3$. The exp map can be expressed in closed form as, with $\mathrm{dexp}_{\mathbf{x}}$ in (A2),

$$ \exp(\hat{\xi}) = \begin{bmatrix} \exp \tilde{\mathbf{x}} & \mathrm{dexp}_{\mathbf{x}}\mathbf{y} \\ \mathbf{0} & 0 \end{bmatrix}. \quad (A6) $$

The vector space $\mathbb{R}^6$ is isomorphic to se (3) by identifying $\hat{\xi} \in se(3)$ with $\xi = (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^6$. Therefore, the exponential map will likewise be written as $\exp(\xi)$ with $\xi \in \mathbb{R}^6$. The vector $\xi$ is the screw coordinate vector that describes a screw motion given by $\exp(\xi) \in SE(3)$. The Lie bracket on se (3) is the screw product, $[\xi_1, \xi_2] = ad_{\xi_1}\xi_2 = (\mathbf{x}_1 \times \mathbf{x}_2, \mathbf{x}_1 \times \mathbf{y}_2 + \mathbf{y}_1 \times \mathbf{x}_2)$, for $\xi_i = (\mathbf{x}_i, \mathbf{y}_i) \in \mathbb{R}^6 \cong se(3)$. The right-trivialized differential of (A6), and its inverse, in terms of canonical coordinates $\xi = (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^6$, are in closed form

$$ \mathrm{dexp}_{\xi} = \begin{bmatrix} \mathrm{dexp}_{\mathbf{x}} & \mathbf{0} \\ (\mathrm{D}_{\mathbf{x}}\mathrm{dexp})(\mathbf{y}) & \mathrm{dexp}_{\mathbf{x}} \end{bmatrix} \quad (A7) $$

$$ \mathrm{dexp}_{\xi}^{-1} = \begin{bmatrix} \mathrm{dexp}_{\mathbf{x}}^{-1} & \mathbf{0} \\ (\mathrm{D}_{\mathbf{x}}\mathrm{dexp}^{-1})(\mathbf{y}) & \mathrm{dexp}_{\mathbf{x}}^{-1} \end{bmatrix} \quad (A8) $$

where $\mathrm{dexp}_{\mathbf{x}}$ and $\mathrm{dexp}_{\mathbf{x}}^{-1}$ are as in (A2) and (A3), and

$$ (\mathrm{D}_{\mathbf{x}}\mathrm{dexp})(\mathbf{y}) = \frac{\beta}{2}\tilde{\mathbf{y}} + \delta(\tilde{\mathbf{x}}\tilde{\mathbf{y}} + \tilde{\mathbf{y}}\tilde{\mathbf{x}}) + \frac{\mathbf{x}^T\mathbf{y}}{\|\mathbf{x}\|^2}(\alpha - \beta)\tilde{\mathbf{x}} + \frac{\mathbf{x}^T\mathbf{y}}{\|\mathbf{x}\|^2}\left(\frac{\beta}{2} - 3\delta\right)\tilde{\mathbf{x}}^2 \quad (A9) $$

$$ (\mathrm{D}_{\mathbf{x}}\mathrm{dexp}^{-1})(\mathbf{y}) = -\frac{1}{2}\tilde{\mathbf{y}} + \frac{1}{\|\mathbf{x}\|^2}(1-\gamma)(\tilde{\mathbf{x}}\tilde{\mathbf{y}} + \tilde{\mathbf{y}}\tilde{\mathbf{x}}) + \frac{\mathbf{x}^T\mathbf{y}}{\|\mathbf{x}\|^4}\left(\frac{1}{\beta} + \gamma - 2\right)\tilde{\mathbf{x}}^2. \quad (A10) $$
The directional derivatives of (A7) and (A8) are

$$ (D_{\xi} \mathrm{dexp})(\eta) = \begin{bmatrix} (D_x \mathrm{dexp})(u) & 0 \\ (D_{\xi} \mathrm{Ddexp})(\eta) & (D_x \mathrm{dexp})(u) \end{bmatrix} \quad (A11) $$

$$ (D_{\xi} \mathrm{dexp}^{-1})(\eta) = \begin{bmatrix} (D_x \mathrm{dexp}^{-1})(u) & 0 \\ (D_{\xi} \mathrm{Ddexp}^{-1})(\eta) & (D_x \mathrm{dexp}^{-1})(u) \end{bmatrix} \quad (A12) $$

where $\eta = [u, v]$. Therein the directional derivatives of matrices $\mathrm{Ddexp}$ and $\mathrm{Ddexp}^{-1}$ possesses the explicitly form

$$ (D_{\xi} \mathrm{Ddexp})(\eta) = \frac{\beta}{2} \tilde{v} + \frac{\alpha - \beta}{\|x\|^2} (x^T u) \tilde{y} + \frac{\alpha - \beta}{\|x\|^2} (x^T y) \tilde{u} + \frac{\beta(2-3\delta)}{\|x\|^2} \left( (x^T u)(\tilde{x}\tilde{y} + \tilde{y}\tilde{x}) + (x^T y)(\tilde{x}\tilde{u} + \tilde{u}\tilde{x}) \right) + \delta (\tilde{x}\tilde{v} + \tilde{v}\tilde{x} + \tilde{y}\tilde{u} + \tilde{u}\tilde{y}) \quad (A13) $$

$$ + \frac{1}{\|x\|^2} \left( (\alpha - \beta)(x^T v + y^T u) - \frac{1}{2}\beta(x^T u)(x^T y) + \frac{1-5\alpha+4\beta}{\|x\|^2}(x^T y)(x^T u) \right) \tilde{x} \\ + \frac{1}{\|x\|^2} \left( \left(\frac{1}{2}\beta - 3\delta\right)(x^T v + y^T u) + \frac{1}{\|x\|^2}(\alpha - \frac{7}{2}\beta + 15\delta)(x^T y)(x^T u) \right) \tilde{x}^2. $$

$$ (D_{\xi} \mathrm{Ddexp}^{-1})(U) = -\frac{1}{2}\tilde{v} + \frac{1}{\|x\|^2} \left( \frac{1}{4}(x^T u)(\tilde{x}\tilde{y} + \tilde{y}\tilde{x}) + (1-\gamma)(\tilde{x}\tilde{v} + \tilde{v}\tilde{x} + \tilde{y}\tilde{u} + \tilde{u}\tilde{y}) \right) \quad (A14) $$

$$ + \frac{1}{\|x\|^4} \left( (x^T u)(\gamma(1+\gamma) - 2)(\tilde{x}\tilde{y} + \tilde{y}\tilde{x}) + (\frac{1}{\beta} + \gamma - 2)(x^T y)(\tilde{x}\tilde{u} + \tilde{u}\tilde{x}) \right. \\ \left. + ((\frac{1}{\beta} + \gamma - 2)(x^T v + y^T u) - \frac{1}{4}(x^T y)(x^T u))\tilde{x}^2 \right) \\ + \frac{1}{\|x\|^6} (x^T y)(x^T u) \left( 8 - \gamma(3+\gamma) - \frac{2}{\beta}(1+\gamma) \right) \tilde{x}^2. $$

The log map on $SE(3)$ is and $SO(3)$

$$ h = \begin{bmatrix} \mathbf{R} & \mathbf{r} \\ \mathbf{0} & 1 \end{bmatrix} \in SE(3) \mapsto \hat{\xi} = \log(h) = \begin{bmatrix} \tilde{x} & \mathrm{dexp}_{\tilde{x}}^{-1} \mathbf{r} \\ \mathbf{0} & 1 \end{bmatrix} \in se(3) \quad (A15) $$

or in vector notation $\xi = [x, \mathrm{dexp}_{x}^{-1} \mathbf{r}] \in \mathbb{R}^6$, where $x \in \mathbb{R}^3$ is obtained with the log on $SO(3)$

$$ \mathbf{R} \in SO(3) \mapsto \tilde{\mathbf{x}} = \log \mathbf{R} = \frac{1}{2 \sin \varphi} (\mathbf{R} - \mathbf{R}^T) \in so(3) \quad (A16) $$

with angle $\varphi = \frac{1}{2} \arccos(\mathrm{tr}(\mathbf{R}) - 1)$.

## B DE CASTELJAU ALGORITHM ON LIE GROUPS

An easy and flexible way to represent 3-dimensional curves in $\mathbb{R}^3$ (and generally on any metric vector space) in terms of control polygons are Bézier curves [6]. The classical De Casteljau algorithm [14, 38] is a computational way to generate the control polygon from a given set of control points [23]. On curved metric spaces, particularly on Riemannian manifolds, including symmetric spaces, and Lie groups, the definition of 'control polygons' must respect the metric structure of the space. Ideally, Bézier curves should be geodesics w.r.t. a certain metric with certain invariance. This complicates the derivation of De Casteljau type algorithms. There is a body of literature on polynomial curves, Bézier curves, and B-splines on (symmetric) Riemannian manifolds and Lie groups [9, 12, 16, 18, 21, 57, 67, 73, 80]. The crucial differences when designing a De Casteljau algorithm on curved spaces is that the subdivision properties and smoothness results of the Bézier curves are not generally inherited from the algorithm on vector spaces [57].

For the special case of compact Lie groups, and particularly for $SO(3)$ and the quaternion group, De Casteljau algorithms where introduced in [55]. This was extended in [25] to spatial motions described by dual quaternions, i.e. the double cover of $SE(3)$, that relies on geodesics w.r.t. a left-invariant metric. The algorithm 9 presented below is an adapted version of [55]. It rests on the assumption that, at the different steps of the algorithm, curves between points are expressed as motion in a 1-parameter subgroup. For the compact rotation group $SO(3)$ these are geodesics, whereas for the non-compact rigid body motion group $SE(3)$ these are screw motions (not geodesics). Possible refinements can be found in [57, 80]. Algorithm 9 constructs a Bézier curve from $n+1$ given points points $h_0, \dots, h_n$ in a Lie group $G$, on which an exp map is defined. The curve interpolates between terminal points $h_0$ and $h_n$, while the remaining points serve as 'control points'. Intermediate points $h_i^k$ are computed recursively.
The Bézier curve of order $n$ is then obtained as $h_0^n(\tau)$. At every $\tau$, its evaluation requires $n$ nested evaluations of the log and exp map, which dictates the computational effort. The algorithm complexity is quadratic in the number $n$ of points. Each $h_i^k$ depends on $\tau$, and symbolic construction of the Bézier curve is generally not possible, even when closed form expressions for exp and log are available. The algorithm is used to evaluate the curve at given $\tau$. The algorithm assumes zero initial and terminal velocities. The case for non-zero boundary velocities is easily incorporated. Details can be found in [16].

**Algorithm 9: De Casteljau Algorithm on Lie Group G**

**1. Input:**

* Samples $h_0, h_1, \dots, h_n \in G$

* Parameter value $\tau \in [0, 1]$

**2. Recursive Computation of 'Control Points':**

FOR $k = 1, \dots, n$ DO

FOR $i = 0, \dots, n-k$ DO

* $\xi_i^k(\tau) := \log((h_i^{k-1})^{-1}h_{i+1}^{k-1})$

* $h_i^k(\tau) := h_i^{k-1}(\tau) \exp(\tau \xi_i^k(\tau))$

END

END

**3. Value of Bézier curve at $\tau \in [0, 1]$:** $h_0^n(\tau)$

**References**

[1] Almaghout, K., Cherubini, A., Klimchik, A.: Robotic co-manipulation of deformable linear objects for large deformation tasks. Robotics and Autonomous Systems (2024).

[2] Almaghout, K., Klimchik, A.: Manipulation planning for cable shape control. Robotics 13(18), 20 (2024)

[3] Altmann, S.: Rotations, Quaternions, and Double Groups. Clarendon Press (1989)

[4] Bartoň, M., Jüttler, B., Wang, W.: Construction of rational curves with rational rotation-minimizing frames via möbius transformations. In: Mathematical Methods for Curves and Surfaces: 7th International Conference, MMCS 2008, Tønsberg, Norway, June 26-July 1, 2008, Revised Selected Papers 7, pp. 15–25. Springer (2010)

[5] Bauchau, O., Epple, A., Heo, S.: Interpolation of finite rotations in flexible multi-body dynamics simulations. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics 222(4), 353–366 (2008)

[6] Bézier, P.: The Mathematical Basis of UNISURF CAD System. Butterworths, London (1986)

[7] Black, C., Till J. and Rucker, D.: Parallel Continuum Robots. IEEE Transactions on Robotics 34(1), 29–47 (2018)

[8] Bloch, A.M.: Nonholonomic mechanics and control. Springer (2003)

[9] Bogfjellmo, G., Modin, K., Verdier, O.: A Numerical Algorithm for $C^2$-Splines on Symmetric Spaces. SIAM Journal on Numerical Analysis 56(4), 2623–2647 (2018)

[10] Bottema, O., Roth, B.: Theoretical Kinematics. North-Holland (1979)

[11] Briot, S., Boyer, F.: A geometrically exact assumed strain modes approach for the geometrico-and kinemato-static modelings of continuum parallel robots. IEEE Transactions on Robotics 39(2), 1527–1543 (2022)

[12] Camarinha, M., Leite, F.S., Crouch, P.: On the geometry of Riemannian cubic polynomials. Differential Geometry and its Applications 15(2), 107–135 (2001)

[13] Carricato, M.: Classification of the three-dimensional persistent poe manifolds of $\mathfrak{se}(3)$. Mechanism and Machine Theory 206, 105926 (2025)

[14] de Casteljau, P.: Enveloppe soleau 40.040. Document no P 2108 (1959)

[15] Condurache, D.: Poisson-Darboux problems's extended in dual Lie algebra. In: AAS/AIAA Astrodynamics Specialist Conference, Stevenson, WA, USA (2017)

[16] Crouch, P., Kun, G., Leite, F.S.: The De Casteljau algorithm on Lie groups and spheres. Journal of Dynamical and Control Systems 5(3), 397–429 (1999)

[17] Crouch, P., Leite, F.S.: Geometry and the dynamic interpolation problem. In: 1991 American Control Conference, pp. 1131–1136. IEEE (1991)

[18] Crouch, P., Leite, F.S.: The dynamic interpolation problem: on Riemannian manifolds, Lie groups, and symmetric spaces. Journal of Dynamical and Control Systems 1(2), 177–202 (1995)
[19] Darboux, G.: Leçons sur la théorie générale des surfaces et les applications géométriques du calcul infinitesimal. Gautiers-Villars, Paris 4 (1887)

[20] Dhullipalla, M.H., Hamrah, R., Warier, R.R., Sanyal, A.K.: Trajectory generation on SE(3) for an underactuated vehicle with pointing direction constraints. In: 2019 American Control Conference (ACC), pp. 1930–1935. IEEE (2019)

[21] Ding, R.: Drawing ruled surfaces using the dual De Boor algorithm. Electronic Notes in Theoretical Computer Science 61, 178–190 (2002)

[22] Duffy, J., Griffis, M., Swinson, M.: The fallacy of modern hybrid control theory for the simultaneous control of force and motion. In: S. Stifter, J. Lenarčič (eds.) Advances in Robot Kinematics, pp. 248–258. Springer Vienna, Vienna (1991)

[23] Farin, G.E.: Curves and surfaces for computer-aided design: A practical guide. Academic Press (1997)

[24] Ge, Q., Ravani, B.: Computer aided geometric design of motion interpolants. Journal of Mechanical Design 116(3), 756–762 (1994)

[25] Ge, Q., Ravani, B.: Geometric Construction of Bézier Motions. Journal of Mechanical Design 116(3), 749–755 (1994)

[26] Hairer, E., Lubich, C., Wanner, G.: Geometric Numerical Integration. Springer (2006)

[27] Han, S., Bauchau, O.A.: On the global interpolation of motion. Computer Methods in Applied Mechanics and Engineering 337, 352–386 (2018)

[28] Hausdorff, F.: Die symbolische Exponentialformel in der Gruppentheorie. Berichte der Königlich-Sächsischen Gesellschaft der Wissenschaften zu Leipzig, Mathematisch-Physische Klasse 58, 19–48 (1906)

[29] Huang, Y., Abu-Dakka, F.J., Silvério, J., Caldwell, D.G.: Toward orientation learning and adaptation in cartesian space. IEEE Transactions on Robotics 37(1), 82–98 (2020)

[30] Huber, P., Perl, R., Rumpf, M.: Smooth interpolation of key frames in a riemannian shell space. Computer Aided Geometric Design 52, 313–328 (2017)

[31] Hussain, I., Malvezzi, M., Gan, D., Iqbal, Z., Seneviratne, L., Prattichizzo, D., Renda, F.: Compliant gripper design, prototyping, and modeling using screw theory formulation. The International Journal of Robotics Research 40(1), 55–71 (2021)

[32] Husty, M.L., Schröcker, H.P.: Algebraic geometry and kinematics. In: I.Z. Emiris, F. Sottile, T. Theobald (eds.) Nonlinear Computational Geometry, The IMA Volumes in Mathematics and its Applications, vol. 151, chap. Algebraic Geometry and Kinematics, pp. 85–107. Springer (2010)

[33] Iserles, A.: Solving linear ordinary differential equations by exponentials of iterated commutators. Numerische Mathematik 45(2), 183–199 (1984)

[34] Iserles, A., Munthe-Kaas, H.Z., Nørsett, S.P., Zanna, A.: Lie-group methods. Acta Numerica 9, 215–365 (2000).

[35] Kang, I., Park, F.: Cubic spline algorithms for orientation interpolation. International journal for numerical methods in engineering 46(1), 45–64 (1999)

[36] Kavan, L., Collins, S., O'Sullivan, C., Zara, J.: Dual quaternions for rigid transformation blending. Trinity College Dublin 5 (2006)

[37] Kim, M.S., Nam, K.W.: Interpolating solid orientations with circular blending quaternion curves. Computer-Aided Design 27(5), 385–398 (1995)

[38] Krautter, J., Parizot, S.: Systeme d'aide la définition eta l'usinage des surfâces de carosserie. Journal de la SIA 44, 581–586 (1971)

[39] Krysl, P., Endres, L.: Explicit newmark/verlet algorithm for time integration of the rotational dynamics of rigid bodies. Int. J. Numer. Methods Eng. 62, 2154–2177 (2005)

[40] Legnani, G., Fassi, I., Tasora, A., Fusai, D.: A practical algorithm for smooth interpolation between different angular positions. Mechanism and Machine Theory 162, 104341 (2021)

[41] Li, Z., Schicho, J., Schröcker, H.P.: Kempe's universality theorem for rational space curves. Foundations of Computational Mathematics 18, 509–536 (2018)

[42] Li, Z., Ye, K.: Rational curves on real classical groups. arXiv preprint arXiv:2408.04453 (2024)

[43] Lovegrove, S., Patron-Perez, A., Sibley, G.: Spline fusion: A continuous-time representation for visual-inertial fusion with application to rolling shutter cameras. In: BMVC, vol. 2, p. 8 (2013)

[44] Magnus, W.: On the exponential solution of differential equations for a linear operator. Communications on pure and applied mathematics 7(4), 649–673 (1954)

[45] Marsden, J.E., Ratiu, T.S.: Introduction to mechanics and symmetry. Physics Today 48(12), 65 (1995)

[46] Marthinsen, A.: Interpolation in lie groups. SIAM Journal on Numerical Analysis 37(1), 269–285 (1999)
[47] Müller, A.: Approximation of finite rigid body motions from velocity fields. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics **90**(6), 514–521 (2010)

[48] Müller, A.: Screw and Lie group theory in multibody dynamics –Motion representation and recursive kinematics of tree-topology systems. Multibody System Dynamics **43**(1), 1–34 (2018).

[49] Müller, A.: Review of the exponential and Cayley map on SE (3) as relevant for Lie group integration of the generalized Poisson equation and flexible multibody systems. Proceedings of the Royal Society A **477**(2253) (2021)

[50] Müller, A.: Fourth-Order Accurate Strain-Parameterized Shape Representation of Beam Elements for Modeling Continuum Robots and Robotic Manipulation of Slender Objects. ASME Journal of Mechanisms and Robotics **17** (2024)

[51] Munthe-Kaas, H.: Runge-Kutta methods on Lie groups. BIT Numerical Mathematics **38**(1), 92–111 (1998)

[52] Munthe-Kaas, H.: High order Runge-Kutta methods on manifolds. Applied Numerical Mathematics **29**(1), 115–127 (1999)

[53] Murray, R., Li, Z., Sastry, S.: A Mathematical Introduction to Robotic Manipulation. CRC Press (1994)

[54] Noakes, L., Heinzinger, G., Paden, B.: Cubic splines on curved spaces. IMA Journal of Mathematical Control and Information **6**(4), 465–473 (1989)

[55] Park, F., Ravani, B.: Bézier curves on Riemannian manifolds and Lie groups with kinematics applications. Journal of Mechanical Design **117**(1), 36–40 (1995)

[56] Park, F.C., Ravani, B.: Smooth invariant interpolation of rotations. ACM Transactions on Graphics (TOG) **16**(3), 277–295 (1997)

[57] Popiel, T.: On parametric smoothness of generalised b-spline curves. Computer aided geometric design **23**(8), 655–668 (2006)

[58] Ravani, R., Meghdari, A.: Spatial rational motions based on rational frenet-serret curves. In: 2004 IEEE International Conference on Systems, Man and Cybernetics (IEEE Cat. No. 04CH37583), vol. 5, pp. 4456–4461. IEEE (2004)

[59] Ren, K., Yuan, T., Liu, J.: Geometrically exact beam finite element with generalized b-spline interpolation on the special euclidean group se (3). Computer Methods in Applied Mechanics and Engineering **441**, 117979 (2025)

[60] Rucker, D.C., Jones, B.A., Webster III, R.J.: A geometrically exact model for externally loaded concentric-tube continuum robots. IEEE transactions on robotics **26**(5), 769–780 (2010)

[61] Sarker, A., Sinha, A., Chakraborty, N.: On screw linear interpolation for point-to-point path planning. In: 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 9480–9487. IEEE (2020)

[62] Selig, J.: Geometric Fundamentals of Robotics. Springer (2005)

[63] Selig, J.: Curves of stationary acceleration in SE (3). IMA Journal of Mathematical Control and Information **24**(1), 95–113 (2007)

[64] Selig, J.: Rational interpolation of rigid-body motions. In: Advances in the Theory of Control, Signals and Systems with Physical Modeling, pp. 213–224. Springer (2010)

[65] Shoemake, K.: Animating rotation with quaternion curves. In: Proceedings of the 12th annual conference on Computer graphics and interactive techniques, pp. 245–254 (1985)

[66] Shoemake, K.: Quaternion calculus and fast animation, computer animation: 3-d motion specification and control. Siggraph (1987)

[67] Sprott, K., Ravani, B.: Kinematic generation of ruled surfaces. Advances in Computational Mathematics **17**, 115–133 (2002)

[68] Tagliavini, A., Bianco, C.G.L.: A smooth orientation planner for trajectories in the cartesian space. IEEE Robotics and Automation Letters **8**(5), 2606–2613 (2023)

[69] Tan, J., Xing, Y., Fan, W., Hong, P.: Smooth orientation interpolation using parametric quintic-polynomial-based quaternion spline curve. Journal of Computational and Applied Mathematics **329**, 256–267 (2018)

[70] Varadarajan, V.S.: Lie groups, Lie algebras, and their representations. Springer Science & Business Media (2013)

[71] Wang, H., Totaro, M., Beccai, L.: Toward perceptive soft robots: Progress and challenges. Advanced Science **5**(9), 1800541 (2018)

[72] Wittenburg, J.: Kinematics. Berlin, Heidelberg: Springer Berlin Heidelberg **10**, 978–3 (2016)

[73] Wu, Y., Carricato, M.: Persistent manifolds of the special euclidean group se (3): A review. Computer Aided Geometric Design **79**, 101872 (2020)

[74] Wu, Y., Müller, A., Carricato, M.: The 2D orientation interpolation problem: a symmetric space approach. In: Advances in Robot Kinematics 2016, pp. 293–302. Springer (2017)

[75] Xiao, Q., Musa, M., Godage, I.S., Su, H., Chen, Y.: Kinematics and stiffness modeling of soft robot with a concentric
backbone. Journal of Mechanisms and Robotics **15**(5), 051011 (2023)

[76] Žefran, M., Kumar, V.: Planning of smooth motions on se (3). In: Proceedings of IEEE International Conference on Robotics and Automation, vol. 1, pp. 121–126. IEEE (1996)

[77] Žefran, M., Kumar, V.: Interpolation schemes for rigid body motions. Computer-Aided Design **30**(3), 179–189 (1998)

[78] Zefran, M., Kumar, V., Croke, C.: Metrics and connections for rigid-body kinematics. International Journal of Robotics Research **18**(2), 243–258 (1999)

[79] Žefran, M., Kumar, V., Croke, C.B.: On the generation of smooth three-dimensional rigid body motions. IEEE Transactions on Robotics and Automation **14**(4), 576–589 (1998)

[80] Zhang, E., Noakes, L.: The cubic de Casteljau construction and Riemannian cubics. Computer Aided Geometric Design **75**, 101789 (2019)