離散メッシュのリーマン幾何を直接扱うロボット学習フレームワーク、DMPとGPとRFMを拡張

Abstract

All the objects composing our world are enclosed within surfaces. Yet, most robot learning and motion generation frameworks treat surfaces as constraints ignoring their intrinsic geometry. This gap is acute for polyhedral meshes—the standard output of CAD and 3D reconstruction—whose discrete geometric structure remains unexploited. In this paper, we propose a unified discrete Riemannian framework that enables robot learning directly on polyhedral surface meshes. Using discrete differential geometry, we define logarithmic and exponential maps, parallel transport, and ambient-space projections that remain well-defined across faces, edges, and vertices. We instantiate the framework in three learning paradigms: (i) Dynamic Movement Primitives (DMPs), an improved exponential-map computation and a fixed-tangent-cone forcing-term encoding with parallel transport yield better cross-surface generalisation and stability over prior mesh-based approaches. (ii) Gaussian Process (GP), a geodesic-based kernel with practical admissibility control, enables regression at arbitrary mesh locations without smoothness assumptions. (iii) Riemannian Flow Matching (RFM), mesh-native operators improve generative quality over spectral baselines while reducing training time. The framework is validated in simulation against state-of-the-art methods and demonstrated on two real-robot scenarios: generalising user-drawn trajectories across different surfaces and planning polishing motions on RGB-D-reconstructed surfaces.

Index Terms:

Learning from Demonstration, Learning on Discrete Manifolds

I Introduction

In real applications, robots have to learn how to physically interact with objects and their environment. Sometimes this is required for manipulation purposes (e.g., grasping, pick-and-place, insertion), and sometimes for manufacturing operations (e.g., wiping, finishing, drawing, grinding, cutting, inspecting). A common trait of all these operations is that they are not executed in free space: in the vast majority of cases, they entail continuous contact with a surface enclosing the target object, where friction effects and geometric irregularities play a significant role. This raises a fundamental question: can a learning process for such tasks be truly effective if the geometric properties of the contact surface are ignored? This issue has been often overlooked in the literature. Most existing Learning from Demonstration (LfD) and motion generation frameworks assume that skills evolve in flat Euclidean spaces and can be directly transferred to surfaces. When task complexity increases, surface constraints are typically recovered a posteriori through projection methods [1], impedance control, or local parameterisation. While effective in structured and repeatable settings, treating surface geometry as an external constraint—rather than an intrinsic component of the learning process—often results in poor generalisation when the underlying structure changes.

In recent years, geometry-aware learning methods have been proposed for some classes of Riemannian manifolds, such as rotations or symmetric positive-definite matrices [2, 3, 4, 5], for which researchers have developed tools and methods to learn robot skills in non-Euclidean spaces. These formulations rely on differential operators—such as exponential and logarithmic maps, and parallel transport—that enable learning and control directly in the tangent space. However, extending these methods to surface-based skills remains an open issue, as surfaces encountered in real applications are rarely available as smooth manifolds. In most cases, surfaces are obtained from CAD models, perception pipelines, or 3D reconstruction tools [6], for which the most natural representation is a 3D polyhedral mesh. Such discrete representations exhibit non-differentiable singularities at edges and vertices and do not admit the smooth structure assumed by classical Riemannian learning frameworks.

The absence of a smooth representation creates a significant gap between geometry-aware robot learning theory and its practical applicability to discrete manifolds. Although surfaces can be studied through the lens of differential geometry, constructing smooth representations for objects with arbitrary topology is generally intractable [7]. At the same time, naïve discretisation of smooth manifold operators on meshes either requires restrictive assumptions on surface smoothness or leads to inconsistencies and ill-defined behaviour at singularities. As a result, existing approaches often rely on local approximations, ad hoc embeddings, or spectral smoothness assumptions that limit robustness and hinder generalisation [8, 9, 10].

Paper Contribution. In this work, we take a more direct route to fill this gap by developing a discrete Riemannian framework for robot learning directly on polyhedral surface meshes. Starting from discrete differential geometry [11], we construct geometrically consistent approximations of key Riemannian operators—including logarithmic and exponential maps, parallel transport, geodesic distance, and ambient-space projections—that remain well-defined across faces, edges, and vertices of a mesh. Unlike prior mesh-based formulations [9], which leave edge and vertex singularities unaddressed and rely on local approximations of the exponential map (valid only for small displacements), our operators are built on the straightest-geodesic principle [12] and are equipped with tangent-cone projections that handle all singular loci. Building on this geometric foundation, we show how three established learning frameworks can be naturally extended to discrete surfaces: (i) geometry-aware Dynamic Movement Primitives (DMPs) with a fixed-tangent-cone forcing-term encoding that replaces the local parameterisations used in [9, 13], (ii) geodesic-distance-based Gaussian Processes (GPs) with controllable kernel admissibility, and (iii) Riemannian Flow Matching (RFM) with mesh-native operators replacing the spectral premetrics of [4]. These formulations yield motion and function representations that generalise across surfaces with different curvature, topology, and discretisation resolution, removing the need for smoothness or spectral assumptions.

To summarize, our main contributions are as follows:

• A discrete Riemannian operator framework on polyhedral meshes—including tangent cones, straightest-geodesic exponential maps, and ambient-space projections—that explicitly handles edge and vertex singularities left unaddressed by existing mesh-based methods [9, 8], enabling consistent geometric computations directly on discrete surface representations like polyhedral meshes.
• A reformulation of geometry-aware DMPs on meshes that encodes forcing term in a fixed tangent cone and reconstructs it via parallel transport during execution, replacing the local, velocity-aligned parameterisation of [9] and eliminating the limit-cycle drift it induces on curved surfaces.
• A geodesic-distance-based GP formulation on meshes that enables kernel regression at arbitrary surface locations—including face interiors—without smoothness or density requirements, overcoming the vertex-restriction and resolution-sensitivity of Laplace–Beltrami-based kernels [10].
• An extension of RFM to discrete polyhedral surfaces using geodesic premetrics, improving generative quality and reducing training time relative to the spectral baselines of [4].

The rest of the paper is organised as follows. Section II analyses the state of the art. After a quick recall of foundation knowledge (Sec. III), we present the unified approach for geometry-aware learning on discrete surfaces in Sec. IV. Section V provides simulation results and comparisons with the state-of-the-art, while Sec. VI provides experiments on a real robot highlighting the strengths of our approach. We conclude in Sec. VII with some limitations of the approach and final remarks.

II Related work

Geometric Learning

Learning from Demonstration (LfD) is a widely adopted paradigm for encoding human expertise and transferring it into robotic systems [14]. Classical LfD approaches—including DMPs [15], Probabilistic MPs [16], Kernelised MPs [17], or Gaussian Mixture Models (GMMs) [18]—have been successfully applied to learn skills defined in Euclidean spaces. However, these methods do not account for the intrinsic geometry of constrained domains, and consequently fail when applied to data evolving on non-linear spaces such as rotations, stiffness matrices, or surfaces. To address this limitation, Lie-group-based formulations and Riemannian geometry have been introduced, with Lie-algebra approaches proving effective for specific spaces [19, 20, 21, 22], and Riemannian formulations offering greater flexibility by imposing fewer structural restrictions on the domain [23, 3, 2, 4, 24]. As a result, geometry-aware learning methods on smooth manifolds have been successfully applied to rotations, symmetric positive-definite matrices, and Grassmannian manifolds, where closed-form expressions of geometric operators are available. However, extending these results to surfaces remains largely an open problem, as real-world surface representations lack the fundamental structure that these methods rely on.

Learning on smooth and discrete manifolds

Several methods have been proposed to parameterise smooth surfaces using spline or subdivision representations [25, 26, 27]. While effective in computer graphics and geometric modelling, these approaches typically provide only local parameterisations or impose strong constraints on surface topology, restricting their applicability to robotic learning tasks involving complex or reconstructed environments. Working directly on polyhedral surface meshes has emerged as a promising alternative. Pantic et al. [8] proposed a mesh-manifold-based motion planning algorithm, framed as a Riemannian motion policy [24], that uses isomorphic mappings to flatten meshes onto planar domains. However, this strategy requires restrictive topological assumptions and implicitly assumes sufficient mesh smoothness, avoiding the singularities arising at edges and vertices rather than addressing them.

The MeshDMP framework [9] represents the most direct prior reference for our approach, as it extends geometry-aware DMPs [23] to discrete surfaces. While MeshDMP removes topological restrictions on the mesh, it suffers from two important limitations. First, it relies on an approximate exponential map that is valid only under bounded face-angle conditions and does not explicitly address vertex singularities. Second, its forcing-term encoding is defined in a local, moving coordinate system, which can lead to instability and degraded generalisation, particularly for rhythmic motion patterns.

Our approach addresses both limitations of MeshDMP and, more broadly, establishes a geometrically consistent foundation for learning on discrete surfaces. By adopting a straightest-geodesic-based [12] exponential map and explicitly handling edge and vertex singularities, we obtain a more principled and consistent geometric formulation. Additionally, by encoding the forcing term in a fixed tangent space and transporting it during execution, we achieve substantially improved stability and transferability across surfaces, without relying on local parameterisations.

Meshes and the Laplace-Beltrami operator

The Laplace-Beltrami operator (LBO) is the generalisation of the Laplacian to smooth Riemannian manifolds. Despite being initially developed for differentiable structures, its definition has been extended to discrete domains, including polyhedral meshes [28, 29]. The LBO is a powerful tool for learning on surfaces, as its spectral decomposition can be used to approximate distances, kernels, and function spaces [30, 31, 32]. As an example, a pioneering work by Borovitskiy et al. [10] exploited the eigenfunction basis of the LBO to construct Matérn kernels on compact manifolds, effectively enabling Gaussian Process Regression (GPR) on meshes.

A key advantage of LBO-based methods is that they offer theoretical guarantees. However, they implicitly assume that the mesh is a sufficiently good approximation of a smooth surface and restrict function evaluation to mesh vertices. As a consequence, accurate GPR often requires dense meshes or additional interpolation steps, increasing computational cost and approximation error. Our approach, instead, relies directly on geodesic distance computation, enabling GPR at arbitrary surface locations—including points in the interior of faces—without any smoothness requirement. It is worth mentioning that geodesic-based kernels may lose positive definiteness beyond a critical lengthscale [33]. We address this limitation by treating it as a practical modelling trade-off and introduce criteria that ensure stable GPR in realistic scenarios.

Riemannian generative models

Continuous Normalising Flows (CNFs) [34] are generative models that represent probability distributions by learning smooth, time-dependent vector fields whose integration defines a continuous transformation from a simple base distribution to a target data distribution. Despite their theoretical appeal and exact likelihood formulation, CNF-based architectures are computationally expensive to train, particularly in high-dimensional settings, and can suffer from numerical instability. Conditional Flow Matching (CFM) [35, 36] addresses these limitations by replacing likelihood maximisation with a regression-based objective, substantially improving both training stability and efficiency.

Particularly relevant for geometric learning is the work by Chen et al. [4], who extend CFM to Riemannian manifolds (RFM) by explicitly incorporating the geometric structure into the flow definition. By defining probability paths on manifolds, their approach trains a neural network to learn vector fields using exact differential operators on analytically defined manifolds (e.g., hyperspheres and tori), or approximate operators for discrete domains such as meshes, constructed via premetrics induced by the Laplace–Beltrami operator.

Our framework extends RFM to polyhedral meshes by introducing geometrically consistent operators tailored to mesh representations, enabling more stable and accurate vector field learning. As demonstrated in the evaluation below, the performance improvement is consistent across the tested baseline methods.

III Background

In this section, we present some background notions needed for a proper understanding of our approach. The notation used throughout the paper is reported in Table I.

III-A Riemannian Manifolds

A Riemannian n-manifold is a pair (\mathbb{M},\mathfrak{g}), where \mathbb{M} is a differentiable manifold and \mathfrak{g}:T_{(\cdot)}\mathbb{M}\times T_{(\cdot)}\mathbb{M}\rightarrow\mathbb{R} is a smoothly varying inner product defined on each tangent space T_{(\cdot)}\mathbb{M} [37]. Based on the definition of the metric \mathfrak{g}, one can construct other meaningful concepts, such as distances, angles, volumes, and curvature.

The choice of \mathfrak{g} plays a critical role and it is application dependent. In this work, we consider the metric to be the one induced by the ambient space of the manifold, i.e.,

\mathfrak{g}:\langle v,w\rangle\triangleq v^{\top}w\quad v,w\in T_{y}\mathbb{M}\subset\mathbb{R}^{2n},

where we define the inner product using the dot product in Euclidean (ambient) space. This is always possible since any smooth n-manifold admits a smooth embedding into \mathbb{R}^{2n} [38], which implies that, for any point p\in\mathbb{M}, the tangent plane T_{p}\mathbb{M} is a n-dimensional linear subspace of \mathbb{R}^{2n}.

Differential geometry is built around the notion that a manifold \mathbb{M} locally resembles its tangent space T_{(\cdot)}\mathbb{M}, which provides a Euclidean approximation at each point. For this reason, specific operators are required to map values from manifolds to tangent spaces, and vice versa. Let p\in\mathbb{M} be a point of interest on the manifold. The logarithmic map \mathrm{Log}_{p}:\mathbb{M}\rightarrow T_{p}\mathbb{M} is the differential operator which enables the expression of any other point on \mathbb{M} in the tangent space T_{p}\mathbb{M} of p. The dual operation of mapping values from the tangent space back to the manifold is achieved through the exponential map \mathrm{Exp}_{p}:T_{p}\mathbb{M}\rightarrow\mathbb{M}. Finally, the parallel transport P_{p\rightarrow q}:T_{p}\mathbb{M}\rightarrow T_{q}\mathbb{M} is the operator that moves vectors from the tangent space of p to the tangent space of another point q\in\mathbb{M}, while preserving their geometric properties (such as length and angle) along the curve connecting p and q on \mathbb{M}. We refer to [2] for the definition, as well as for a good visual and intuitive representation, of these operators for popular manifolds in robotics.

III-B Polyhedral Meshes

A polyhedral mesh is a discrete representation of a surface formed by a finite collection of planar polygonal faces, typically triangles, joined along shared edges and vertices. This representation is widely used in Computer-Aided Design (CAD) and computer vision, and is supported by a wide range of data structures and algorithms available off-the-shelf in the form of libraries [39].

In this work, we focus on manifold surface meshes \mathcal{M}, which can be viewed as piecewise-linear, two-dimensional manifold surfaces embedded in \mathbb{R}^{3}. In contrast to the smooth surfaces of differential geometry, polyhedral meshes are flat on each face and exhibit discontinuities in the normal field at edges and vertices, where tangent planes are not uniquely defined. For this reason, polyhedral meshes are commonly referred to as discrete surfaces.

A mesh \mathcal{M} is described as a graph whose nodes are the vertices \nu_{i}\in\mathcal{V}, connected by edges, and whose faces \mathcal{T}_{i}\in\mathcal{F} are the planar polygons bounded by those edges; without loss of generality, we consider faces to be triangles. Relevant to our work is the definition of the total angle \theta(p) for points p\in\mathcal{M}. Following [12], the total angle for points which are either in the interior or on the edge of a triangular face is \theta(p)\triangleq 2\pi, since the surface is locally isometric to the Euclidean space at those points; at a vertex \nu_{i}, we define

\theta(\nu_{i})\triangleq\sum_{\mathcal{T}_{j}\in\mathcal{I}(\nu_{i})}\theta_{\mathcal{T}_{j}}(\nu_{i}),

where \mathcal{I}(\nu_{i}) is the set of all faces adjacent to vertex \nu_{i}, and \theta_{\mathcal{T}_{j}}(\nu_{i}) is the angle subtended by face \mathcal{T}_{j} at vertex \nu_{i}. Based on the total angle, points are classified as spherical (\theta(\nu_{i})<2\pi), Euclidean (\theta(\nu_{i})=2\pi), or hyperbolic (\theta(\nu_{i})>2\pi), as shown in Fig. 1.

III-C Geometry-Aware Dynamic Movement Primitives

Dynamic Movement Primitives (DMPs) are an established framework to learn a generalisable representation of sensorimotor skills from demonstrations, by encoding expert demonstrations into a dynamical system [40, 15]. Initially developed to operate in Euclidean space, DMPs have sparked interest throughout the robotics community, which has proposed their extension to different non-Euclidean manifolds.

Among different available formulations, we recall the discrete geometry-aware DMP [23] which uses the Riemannian operators on smooth manifolds:

\displaystyle\tau\nabla_{z}z \\ \displaystyle=\alpha(\beta\mathrm{Log}_{p}(g)-z)+f(x), \\ \displaystyle\tau\dot{p} \\ \displaystyle=z, \\ \displaystyle\tau\dot{x} \\ \displaystyle=-\zeta x.

Here, is the position state (point on the manifold), g\in M is the dynamics attractor, z\in T_{p}M is the scaled velocity state, x\in\mathbb{R} is a time encoding variable, and f:\mathbb{R}\rightarrow T_{p}M is the forcing function which is learned from a demonstration, and is often modelled as a radial basis functions network.

III-D Gaussian Processes on Manifolds

GPs are powerful and flexible probabilistic models that define, in a principled manner, distributions over functions conditioned on observed data [41]. Let h:M\rightarrow\mathbb{R} be a real-valued random function defined over a manifold M1. The function h is said to follow a GP with mean function \mu and covariance function k, i.e.,

h\sim\mathcal{GP}(\mu,k),

if, for any finite set of m input points P=\{p_{i}\in M\}_{i=1}^{m}, the corresponding function values H=h(P)\in\mathbb{R}^{m} are jointly Gaussian distributed with mean \mu_{P}=\mu(P)\in\mathbb{R}^{m} and covariance matrix K_{PP}=k(P,P)\in\mathbb{R}^{m\times m}. The key component of a GP is the kernel function k:M\times M\rightarrow\mathbb{R}, which encodes prior information and assumptions about the smoothness and structure of the underlying function.

Gaussian process regression (GPR) is the technique that exploits GPs to perform principled probabilistic predictions of the function value outside the training dataset. Given some observed input-output pairs (p_{i},y_{i}), with y_{i}=h(p_{i})+\varepsilon_{i} and \varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2}), and a prior h\sim\mathcal{GP}(0,k), the posterior distribution h\mid Y is another GP with the following mean and covariance:

\displaystyle\mathbb{E}\left\{h\mid Y\right\} \\ \displaystyle=K_{(\cdot)P}\big(K_{PP}+\sigma^{2}I\big)^{-1}Y, \\ \displaystyle\mathbb{C}\textrm{ov}\left\{h\mid Y\right\} \\ \displaystyle=K_{(\cdot,\cdot)}-K_{(\cdot)P}\big(K_{PP}+\sigma^{2}I\big)^{-1}K_{P(\cdot)},

with Y=(y_{i}) and (\cdot) denoting arbitrary test locations.

III-E Riemannian Flow Matching

RFM [4] is a generative modelling framework that extends Continuous Normalising Flows (CNFs) to data supported on Riemannian manifolds in a geometrically consistent manner. CNFs learn a smooth, invertible transformation that maps a simple reference distribution into a complex target distribution; RFM generalises this idea to curved spaces by ensuring that the transformation respects the geometry of the manifold. Let \pi:M\rightarrow\mathbb{R}^{+} be a probability density function on a manifold, and let \mathcal{P} be the space of probability densities. RFM seeks to find a probability path \pi_{t}:[0,1]\rightarrow\mathcal{P} that continuously deforms a simple reference distribution \pi_{\textrm{ref}}=\pi_{0} into a target data distribution \pi_{\textrm{data}}=\pi_{1}. Such a path is induced by the flow \psi_{t}:M\rightarrow M of a time-dependent vector field u_{t}:\mathbb{R}\times M\rightarrow T_{(\cdot)}M, defined as the solution of

\frac{\mathrm{d}}{\mathrm{d}t}\psi_{t}(p)=u_{t}\big(\psi_{t}(p)\big),\quad\psi_{0}(p)=p.

Intuitively, \psi_{t} moves particles on the manifold along trajectories dictated by the vector field u_{t}, progressively transporting mass from \pi_{\textrm{ref}} to \pi_{\textrm{data}}. The RFM objective is to learn a parametric vector field v_{t}(p,\Theta)\in T\mathbb{M} —generally using Neural Networks (NNs)— which best approximates u_{t} by minimising

\mathcal{L}_{\textrm{RFM}}(\Theta)=\mathbb{E}_{t\sim\mathcal{U}(0,1),p\sim\pi_{t}}\left\|v_{t}(p,\Theta)-u_{t}(p)\right\|^{2}_{\mathfrak{g}},

where the norm \|\cdot\|_{\mathfrak{g}} is induced by the Riemannian metric, ensuring that the discrepancy between vector fields is measured in a geometrically meaningful way. Once v_{t}(\cdot,\cdot) is learned, samples from \pi_{\textrm{data}} are generated by integrating the flow starting from p\sim\pi_{\textrm{ref}}. The main challenge lies in the construction of the supervision signal u_{t}. In FM, a premetric D:M\times M\rightarrow\mathbb{R} is a bivariate function on the manifold which is positive definite (D(p,q)\geq 0\ \forall\,p,q\in M and D(p,q)=0\Leftrightarrow p=q) and non-degenerate (\nabla_{p}D(p,q)\neq 0\Leftrightarrow p\neq q); intuitively, D generalises the notion of distance without requiring the full structure of a metric. Given a monotonically decreasing scheduler \eta:[0,1]\rightarrow\mathbb{R} with \eta(0)=1 and \eta(1)=0, the FM objective defines a flow \psi_{t} that shrinks the premetric D(\cdot,p_{1}) according to the law

D\big(\psi_{t}(p_{0}|p_{1}),p_{1}\big)=\eta(t)D(p_{0},p_{1}),

meaning that particles are steered monotonically closer to their target p_{1} as t progresses. Chen et al. [4] show that such a flow can be obtained by defining the conditional vector field u_{t} as

u_{t}(p|p_{1})=\frac{\mathrm{d}\log\eta(t)}{\mathrm{d}t}\,D(p,p_{1})\frac{\nabla D(p,p_{1})}{\|\nabla D(p,p_{1})\|^{2}_{\mathfrak{g}}}.

IV Methodology

IV-A Discrete differential operators

In this work, we aim at extending methods defined for Riemannian manifolds to polyhedral meshes. For this reason, it is necessary to investigate how to properly define the operators described in Sec. III-A, namely logarithmic and exponential maps, and parallel transport, in the space of meshes \mathcal{M}.

IV-A1 Tangent cone

On smooth manifolds, it is always possible to define a locally linear approximation of \mathbb{M} at a point p, the tangent space T_{p}\mathbb{M}. Due to the discrete nature of meshes \mathcal{M}, this concept is ill-defined at edges and vertices. Still, since meshes are Lipschitz manifolds, we can rely on the natural generalisation of tangent spaces: tangent cones. Following [42], let p_{i} be a sequence of points on \mathcal{M} converging to p, and t_{i} a positive sequence decreasing to 0; then the tangent cone T_{p}\mathcal{M} consists of all vectors v expressible as:

v=\lim_{i\rightarrow\infty}\frac{p_{i}-p}{t_{i}}.

From this definition, we may identify three cases for tangent cones, as depicted in Fig. 2. If the point p lies in the interior of a face \mathcal{T}_{i}, then T_{p}\mathcal{M} is isomorphic to the 2-dimensional plane containing \mathcal{T}_{i}. If p lies on an edge, the tangent cone consists of the union of the two half-planes defined by the incident faces on the shared edge. Finally, the tangent cone at a vertex v_{i} is formed by all directions spanned by the faces incident to v_{i}, resulting in a generally non-flat, piecewise-linear cone.

IV-A2 Logarithmic map

The logarithmic map

\displaystyle\text{Log}_{p}:\,\mathcal{M} \\ \displaystyle\rightarrow T_{p}\mathcal{M} \\ \displaystyle q \\ \displaystyle\mapsto v

is the operator that maps a point q on the mesh into a vector v in the tangent cone at p. The vector v represents the initial direction of the shortest path (geodesic) \gamma from p to q, with a magnitude \|v\| equal to the geodesic distance d_{\mathcal{M}}(p,q).

As in the smooth case, this operation requires finding the shortest path \gamma, parameterised by arc-length, joining points p and q. While this problem can be challenging for arbitrary smooth manifolds, robust algorithmic solutions for the shortest path problem exist for polyhedral meshes. In this respect, the continuous Dijkstra algorithm [43] serves as foundation work to compute discrete geodesics \gamma. This method generalises the classical graph-based Dijkstra algorithm to continuous domains by propagating a wavefront over the faces of a triangulated surface, maintaining exact Euclidean distances while respecting surface constraints. In this work, we exploit enhanced versions of this algorithm [44, 45] for more efficient implementations. These methods are based on constructing offline accelerated query structures that enable the subsequent real–time generation of shortest paths in fractions of a millisecond.

The shortest path \gamma on \mathcal{M} connecting p and q is a polyline which can be easily parameterised by arc-length in the domain I=[0,\ell], with \ell being the path length. In view of this definition, the vector representation v of the point q in T_{p}\mathcal{M} given by the logarithmic map is defined as

v=\mathrm{Log}_{p}(q)\triangleq\ell\,\dot{\gamma}(0).

Figure 3 reports a visual representation of this operator.

IV-A3 Exponential map

The exponential map

\displaystyle\mathrm{Exp}_{p}:\,T_{p}\mathcal{M} \\ \displaystyle\rightarrow\mathcal{M} \\ \displaystyle v \\ \displaystyle\mapsto q

serves as the inverse operation to the logarithmic map; it maps a vector v, defined in the tangent space (or tangent cone) of a point p, to a point q on the mesh \mathcal{M}.

The exponential map is related to the solution of the initial value problem (IVP) for the geodesic path \gamma. Parameterised by arc-length, \gamma starts at \gamma(0)=p with initial velocity \dot{\gamma}(0)=v/|v| and has an overall length \ell=|v|. In smooth manifolds, the solution to this initial value problem is always (locally) unique [46], and the geodesic path can be computed via the integration of an Ordinary Differential Equation (ODE) that depends on the intrinsic curvature (and its derivative) of the manifold [47].

In the discrete setting of meshes, the ODE integration approach is ill-defined, as derivatives cannot be evaluated at edges and vertices, which constitute the non-differentiable points of the manifold. A possible strategy to address edge singularities is to compute the planar unfolding of adjacent faces to propagate the vector [43]. Despite this improvement, the solution to the IVP remains ambiguous at vertices, as no unique shortest path may pass through spherical vertices, and an infinite number of shortest paths may pass through hyperbolic points. A unified approach to addressing both vertex and edge singularities was proposed by Polthier and Schmies [12]. In particular, they noted that smooth geodesics are shortest paths, in the sense that they minimise length, but are simultaneously straightest paths, in the sense that they have zero geodesic curvature when measured relative to the manifold itself. Whilst the former property is insufficient to uniquely determine the initial value problem at mesh vertices [43], the latter consistently yields a unique solution.

Let \gamma\subset\mathcal{M} be a curve on the mesh; given a point g\in\gamma, it is possible to compute the left curve angle \theta_{l} (respectively, right angle \theta_{r}) at the point as the sum of the angles, measured on the mesh, covered clockwise (respectively, counter-clockwise) from the vectors with directions -\dot{\gamma}(g) and \dot{\gamma}(g), as shown in Fig. 4. By [12, Def. 6], we say that \gamma is a straightest geodesic if, for any point g in the curve, the left and right curve angles \theta_{l} and \theta_{r} are equal, and evaluate to half the total angle (2).

Exploiting this definition, we define a numerical procedure—outlined in Algorithm 1—to compute the exponential map (14) that is always well-defined, even when v passes through vertices of \mathcal{M}. In a iterative fashion, the algorithm displaces p along the direction v and uses the straightest geodesic definition to propagate motion over non-differentiable points for the manifold. Upon initialisation (line 1), we compute the tentative displaced point \hat{p}_{k+1} (line 3). If such a point lies in the interior of the face \mathcal{T}_{k} containing p_{k} and v_{k} (line 4–5), then we completed the straightest geodesic propagation and \hat{p}_{k+1} represents the final exponential mapping q (line 6–7); this case is visually depicted in Fig. 5(a). Otherwise, Fig. 5(b), when \hat{p}_{k+1} falls outside \mathcal{T}_{k}, we apply straightest geodesic propagation at the singular point; let p_{k+1} be the point in which v_{k} intersects its containing face (line 9), we compute its total angle \theta (line 10) and the remainder vector w (line 11) which accounts for the motion between p_{k} and p_{k+1}; by rotating -w intrinsically on the mesh by \theta/2 (line 12), we construct the vector v_{k+1} generating the straightest path at p_{k+1}. By iterating this procedure, we ensure that v_{k} constantly decreases in length, finally reaching the terminal projection point q.

IV-A4 Parallel transport

The parallel transport

\displaystyle P_{p\rightarrow q}:\,T_{p}\mathcal{M} \\ \displaystyle\rightarrow T_{q}\mathcal{M} \\ \displaystyle v \\ \displaystyle\mapsto w

is the operator that enables the expression of a vector v, defined in the tangent cone of a point p, into the different tangent cone placed at q. More in detail, P_{p\rightarrow q} defines an isometric transformation that moves vector v into w with zero covariant derivative along the path connecting p and q.

In our case, considering that the mesh \mathcal{M} is embedded in the Euclidean space \mathbb{R}^{3}, it follows that the isometric transformation is a rotation, that we can encode as a matrix R\in\mathbb{R}^{3\times 3}. To construct such a rotation, we first need to solve the boundary value problem of computing the geodesic path \gamma connecting p and q, as discussed in Sec. IV-A2. Since the velocity vector field is parallel to itself when moving along a geodesic, the parallel transport transformation (15) can be evaluated as

\displaystyle R(\ell) \\ \displaystyle\triangleq\begin{bmatrix}\dot{\gamma}(\ell)&n(\ell)\times\dot{\gamma}(\ell)&n(\ell)\end{bmatrix}, \\ \displaystyle w \\ \displaystyle\triangleq R(\ell)R^{\top}(0)\,v\in T_{q}\mathcal{M},

where n(\cdot):I\rightarrow\mathbb{R}^{3} is function that assigns each point on the geodesic \gamma its surface normal vector; this operator is visually represented in Fig. 6.

IV-A5 Ambient space and projections

Most learning methods are based on learning some vector fields \psi:\mathcal{M}\rightarrow T_{(\cdot)}\mathcal{M} defined in the manifold itself. The main problem of this formulation is that it hardly integrates with existing, state of the art implementations which are defined in Euclidean spaces.

For this reason, being able to work in the ambient space \mathbb{R}^{3} and constructing an ”augmented” vector field \psi^{\circ}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3} approximating \psi on \mathcal{M} is often a simpler approach. Based on this idea, it is necessary to define operations to project points and vectors back and forth between the mesh and its embedding space, as depicted in Fig. 7.

While transforming a point p from the mesh to p^{\circ} on Euclidean space is trivial, the inverse operation is not. The simplest approach is, in this case, to use the closest point function \Pi_{\mathcal{M}} to project point back to the discrete manifold:

\Pi_{\mathcal{M}}(p^{\circ})\triangleq\arg\min_{p\in\mathcal{M}}\|p-p^{\circ}\|.

In practice, this projection can be efficiently computed using suitable data structure—such as Axis-Aligned Bounding-Boxes [48]—that accelerate closest-point queries.

In the same way, for vectors in the tangent cone T_{p}\mathcal{M}, which are trivially embedded in \mathbb{R}^{3}, we define \textrm{Proj}_{p} as the projection function that maps vector in \mathbb{R}^{3} back to the tangent cone T_{p}\mathcal{M} at the point p:

\textrm{Proj}_{p}(v^{\circ})\triangleq\arg\min_{v\in T_{p}\mathcal{M}}\|v-v^{\circ}\|.

If p lies in the interior of a face, then T_{p}\mathcal{M} becomes a plane and \textrm{Proj}_{p}(v^{\circ}) simplifies to the normal projection to the plane itself. In all other cases, we rely on the numerical procedure described in Algorithm 2 to perform the projection: first, we project the vector v normally onto all faces (lines 2–8) and onto directed edges (lines 9–14) incident to the point p, and remove all projections outside the tangent cone T_{p}\mathcal{M}; we then select (line 15) as projection the candidate which maximises the dot-product w.r.t. the initial vector v.

IV-B G-DMPs on Meshes

From the above formulation of the operators on the mesh manifold, it is straightforward to generalise the geometry-aware DMP formulation in (3) to the space of discrete meshes. The core elements that enables the learning of the actual demonstration in DMP is the forcing term f, a function that shall map the phase variable in \mathbb{R} into the tangent cone T_{p}M of the current position state p. A standard way of learning such forcing term is by considering the tangent bundle as an inclusion of the ambient space of the manifold itself; in this way, the DMP effectively learns a forcing term living in a bigger space dimension, e.g., \mathbb{R}^{3} in the case of surfaces, and then either a) they project f in T_{p}M before integrating the systems dynamic, or b) they integrate the unconstrained dynamical system and perform projection of the position state y once it starts deviating too much from the manifold itself.

We believe that such approach yields poor generalisation capability: considering a policy learned in the trivial surface given by the xy plane, then any attempt in integrating the same forcing term on the xz plane would not produce meaningful results. To overcome this issue, other works [13, 9] defined the forcing term f to be locally parameterised by the current state p; while this solution overcomes the problem of having transferable motion policies, the locality of the parameterisation may fail in capturing global behaviour of the motion policy itself.

In this work, we propose to learn the forcing term f in the tangent cone of a fixed point, and then use parallel transport when executing the motion to retrieve a value properly defined in T_{p}\mathcal{M}; for this scope, we select as fixed point the goal configuration g. Given a demonstration \mathcal{D}=\left\{p_{k},\dot{p}_{k},\ddot{p}_{k}\right\}_{k=1}^{D}, with p_{k}\in\mathcal{M} and \dot{p}_{k},\ddot{p}_{k}\in T_{p}\mathcal{M}, firstly we compute the local forcing term f_{p,\cdot} by inverting (3), and then use parallel transport to compute the forcing term f_{g,\cdot} at the motions goal:

\displaystyle f_{p,k} \\ \displaystyle=\tau^{2}\ddot{p}_{k}-\alpha\big(\beta\mathrm{Log}_{p_{k}}(g)-\tau\dot{p}_{k}\big), \\ \displaystyle f_{g,k} \\ \displaystyle=P_{p_{k}\rightarrow g}\left(f_{p,k}\right).

In this way, the function f_{g,k} lies in T_{g}\mathcal{M} for all D samples; upon the definition of an isomorphism from the bi-dimensional subspace in \mathbb{R}^{3} given by T_{g}\mathcal{M} to \mathbb{R}^{2}, the forcing term can be encoded in a minimal representation as the linear combination of N_{b} non-linear basis functions.

To generalise the motion policy to other final points or different surfaces, it is sufficient to select a new g^{\prime}\in\mathcal{M}^{\prime} in a potentially different mesh \mathcal{M}^{\prime}, and choose a suitable parameterisation to reconstruct a forcing term f\in T_{g^{\prime}}\mathcal{M}^{\prime} from the encoded 2D representation. While integrating, f(x) is first evaluated at T_{g^{\prime}}\mathcal{M}^{\prime} and then parallel transported to the current position p tangent cone T_{p}\mathcal{M}^{\prime}, effectively yielding the following manifold-ODE:

\displaystyle\tau\nabla_{z}z \\ \displaystyle=\alpha(\beta\mathrm{Log}_{p}(g)-z)+P_{g\rightarrow p}\big(f(x)\big), \\ \displaystyle\tau\dot{p} \\ \displaystyle=z, \\ \displaystyle\tau\dot{x} \\ \displaystyle=-\zeta x.

While (22) can be used to encode and generalise discrete motions, we can use the rhythmic formulation of DMPs to enable learning and execution of periodic motion:

\displaystyle\nabla_{z}z \\ \displaystyle=\Omega\Big(\alpha(\beta\mathrm{Log}_{p}(g)-z)+P_{g\rightarrow p}\big(f(\phi)\big)\Big), \\ \displaystyle\dot{p} \\ \displaystyle=\Omega z, \\ \displaystyle\dot{\phi} \\ \displaystyle=2\pi\Omega,

where in this case g becomes the centre of the motion and \Omega is a time scaling factor associated to the frequency of execution of the motion policy.

IV-C Geodesic-based Gaussian Processes

In this work, we are interested in evaluating GPs (4) for which the input space \mathcal{X} is the mesh \mathcal{M}. When working on manifolds, the main challenge of using GPs is the definition of a valid kernel function k that must be positive definite for any choice of points in the manifold itself. In Euclidean domain, the squared exponential kernel

k(x,x^{\prime})=\sigma^{2}\exp\left(-\frac{\|x-x^{\prime}\|^{2}}{2\kappa^{2}}\right)

is one of the most popular choices to induce infinitely mean-square differentiable GPs, with the variance \sigma^{2} and length scale \kappa as hyperparameters.

In this work, we propose to use a kernel in which the norm \|\cdot\| in the Euclidean space is simply replaced with the geodesic distance d_{g} on the manifold itself:

k_{d}(x,x^{\prime})=\sigma^{2}\exp\left(-\frac{d_{g}^{2}(x,x^{\prime})}{2\kappa^{2}}\right).

While [33, Theorem. 2] proves that such a kernel may yield negative-definite covariance matrices for some \kappa>\kappa_{th}, we believe that in most practical application the threshold \kappa_{th} is big enough to guarantee a correct GPR.

For this reason, we designed Algorithm 3 to estimate the value of the maximum admissible length scale value given the set of training data for the GP. We first compute the m\times m distance matrix D between the a set of input points (lines 1–2); then, we follow a bisection based approach to find the critical length scale \kappa_{th}: let \kappa_{min} (respectively, \kappa_{max}) be a length scale for which the kernel is positive (respectively, negative) definite, we first initialise \kappa_{min}=\kappa_{0} and \kappa_{max}=2\kappa_{0} and iterate (lines 3–9) until we find the pair of length scales satisfying our assumptions; we proceed (lines 10–17) using a standard bisection algorithm when we iteratively evaluate and update the search range based on the positive definiteness of the kernel at the length scale \kappa_{mid} in the middle of \kappa_{min} and \kappa_{max}, until numerical convergence. Algorithm 3 can be used to compute the critical length scale \kappa_{th} when the set P of test point is known; if such set is not available upfront, it is also possible to run a Montecarlo identification of \kappa_{th} by randomly sampling points on \mathcal{M}, computing the critical length scale with the algorithm, and choosing the minimum \kappa_{th} among all tests.

IV-D Riemannian Flow Matching on Meshes

To successfully implement RFM on meshes, it is necessary to choose a valid premetric D for the space of discrete surfaces. For general geometries, original work [4] proposes the use of spectral-based distances—which can be computed by the Laplace-Beltrami operator—and use a ODE solver on that metric to compute the reference vector field of the training objective (8).

In this work, we propose to directly use the shortest path on \mathcal{M} to measure distance between samples, effectively choosing the geodesic distance d_{g} as valid premetric. Let p_{0}\backsim\pi_{\textrm{ref}} and p_{1}\backsim\pi_{\textrm{data}} two samples on the mesh drawn respectively from the reference and data distributions, we select as corresponding flow \psi_{t} the shortest path \gamma:[0,1]\rightarrow\mathcal{M} where \gamma(0)=p_{0} and \gamma(1)=p_{1}; once \gamma is constructed, sampling points p_{t} along the flow for times t\in\mathcal{U}(0,1), as well as the corresponding tangential velocity vector u_{t}\propto\dot{\gamma}(t), becomes trivial. Algorithm 4 shows how to train RFM on meshes, using the linear scheduler \eta(t)=1-t as in [4].

To avoid restrictions on the topology of the mesh and learn a parameterisation-free vector field, the NN describing v_{t} directly works in the ambient space \mathbb{R}^{3} of the mesh. Point from the reference distribution \pi_{\textrm{ref}} are uniformly sampled by first selecting a mesh face with probability proportional to its area, and then drawing a point uniformly within the interior of that face.

To drawn a sample distributed as \pi_{\textrm{data}} given the trained vector field v_{t}, we first extract a sample p_{0} from the reference distribution and perform a mesh-based integration technique over the field by consistently projecting the NN outcome on the mesh tangent cone using (19), followed by the geometrically-consistent exponential map operation (14), as shown in Algorithm 5.

V Simulation Results

In this section, we evaluate and compare the proposed methodology against state-of-the-art approaches. Section V-A compares our DMP with the MeshDMP framework [9], Sec V-B compares our geodesic-based squared exponential kernel with the one computed from the Laplace-Beltrami operator introduced in [10], while Sec. V-C compares our Mesh Flow Matching against the RFM approach in [4].

Our method has been implemented fully in C++ and it relies on the CGAL library for polyhedral mesh operations [49, 50]. Experiments in Sec. V-A and V-B have been carried out on a commercial laptop equipped with an Intel Core i7 Ultra 155H, 24GB of RAM, and a NVidia RTX 4060 graphics card with 8GB of VRAM. For the experiments on Flow Matching (Sec. V-C), we used a workstation equipped with an AMD Ryzen 9900X CPU with with 12 cores, 96GB of RAM, and a NVidia RTX 5090 graphics card with 32GB of VRAM.

For testing, we used several meshes, whose main properties are summarised in Table II.

V-A Comparison against MeshDMP

We compare against the MeshDMP [9] framework on the representative task of learning an 8-shaped trajectory demonstrated in the trivial, flat surface, and generalise the policy to the Stanford Bunny mesh. Both approaches utilise same DMP parameters (\alpha=48, \beta=\alpha/4, \gamma=5, and 20 basis) and are initialised at the same initial state, allowing the direct behaviour comparison between MeshDMP and the rhythmic DMP (23).

As depicted in Fig. 8, both methods quickly converge to the limit cycle induced by the 8-shaped trajectory in the initial phase of execution. Differences, however, arise when we let system integrate their dynamics over a longer time horizon. Particularly, we observe that DMP (23) remains faithful to the limit cycle, while MeshDMP exhibits an undesired rotational behavior around the fixed motions centre. Such a difference can be attributed on how the two methods encode the forcing term f used during DMP integration. MeshDMP encodes this term using a local parameterisation that fails in capturing the dynamics change induced by the mesh curvature. Our method, instead, overcomes this limitation by encoding f in a fixed coordinate system and use a geometrically-consistent parallel transport operation to ensure that the forcing term is correctly applied at the actual DMP position state.

Additionally, other known limitations of MeshDMP may contribute to these behaviours. Remarkably, Dalle Vedove et al. relied on a local approximation of the exponential map based on surface normal projection, which is only accurate for sufficiently small vectors and restricts its applicability to surfaces whose relative angle between edge-adjacent faces is less then 90^{\circ}. The local approximation of MeshDMPs exponential map is shown in Fig. 9, where we compare the two integration schemes applied to a long vector, showing how exponential map (14) can faithfully reconstruct a geodesic on the Mesh.

V-B Gaussian Processes

In this section, we evaluate the performance of the proposed geodesic-based kernel for GPR on discrete surfaces. We compare our approach against the state-of-the-art work by Borovitskiy et al. [10], where authors used the spectral properties of the Laplace-Beltrami operator to define kernels on compact manifolds, including the discrete space of meshes. Starting by the consideration that, in Euclidean space, GPs are the solution of a stochastic differential equations that depends on the Laplacian \Delta [53], they generalise such formulation to Riemman manifolds by replacing the Laplacian with the Laplace-Beltrami operator \Delta_{g}; by computing the orthonormal basis given by the eigenfunction-eigenvalue pairs \{(f_{n},\lambda_{n})\} such that

-\nabla_{g}f=\sum_{n=0}^{\infty}\lambda_{n}\langle f,f_{n}\rangle f_{n},

authors define the squared exponential kernel to be

k_{\Delta}(x,x^{\prime})=\frac{\sigma^{2}}{C}\sum_{n=0}^{\infty}e^{-\frac{\kappa^{2}}{2}\lambda_{n}}f_{n}(x)f_{n}(x^{\prime}).

Here, \kappa and \sigma are respectively the length scale and variance hyperparameters as in (25), and C is a normalising constant; for computational purposes, the evaluation of the kernel (27) is truncated.

While the kernel (27) offers the advantage of being positive-definite for any value of \kappa, it suffers from at least two major limitations: a) it implicitly assumes the mesh to be sufficiently smooth, effectively considering it as an approximation of a surface; b) the Laplace-Beltrami operator \Delta_{g} can only be evaluated at the vertices of such a mesh. It follows that predicting function values at arbitrary points, specifically ones located at the interior of the faces, either requires additional interpolation steps, thus introducing approximation errors, or upsampling the mesh to let the test point become a new vertex for the mesh, effectively increasing the model complexity. The geodesic-based kernel (25), instead, automatically solves this problem by enabling the computation of geodesics at arbitrary points.

As an illustrative example, we consider the problem of performing GPR on a non-smooth, trivial surface composed by 3 vertex-adjacent faces of a cube, shown in Fig. 10. A minimal representation of this mesh consists of 7 vertices and 3 faces2; however, to enable a valid application of the Laplace-Beltrami-based kernel (27), we rely on midpoint subdivision to upsample the mesh up to 817 vertices and 1\,536 faces. We additionally test GPR on a further upsampled mesh with 12\,481 vertices and 24\,576 faces. To evaluate k_{\Delta}, we used the GeometricKernels package [54], and use the first 1\,000 eigenfunction as basis. By constructing a ground-truth function on the mesh, we compute a noise-free GP prediction (5) using both kernels, and evaluate regression performance using the mean absolute error (MAE) metric. We run 100 experiments, each randomly selecting 70 vertices from the upsampled mesh as training data, enabling a statistical characterisation of performance across different length scale values. While the two methods exhibit comparable performance at small length scales, we observe that the Laplace-Beltrami-based kernel becomes unstable for larger values of \kappa, with the MAE in some cases degrading to unacceptable levels. In contrast, the geodesic-based kernel remains well behaved, provided that the critical length \kappa_{th} is not reached. Moreover, as shown in Fig 10(d), further increasing the polygonal density of the mesh, which is used to discretise the Laplace-Beltrami operator on a non-smooth surface, leads to degraded performance. This behaviour is a consequence of the assumption of a smooth underlying manifold implicitly made by the Laplace-Beltrami operator: increasing the mesh resolution in the presence of geometric discontinuities, amplifies the discretisation errors, ultimately resulting in a breakdown of the GPR model.

A similar experiment, Fig. 11, has been ran on the Stanford Dragon mesh with 9\,983 faces and 4\,898 vertices. Also in this case, our proposed kernel performs slightly better than the Laplace-Beltrami based solution in the whole spectrum of length scales below the critical value of k_{th}=0.023, which is assessed through the Algorithm 3.

Computation performance

To fully compare the methods, we also evaluate their computational performance both in terms of time required to construct the necessary data structures and to compute the covariance matrices. The summarising results are reported in Table III.

Within our method, it is necessary to construct a data structure for each set of input point for the kernel; upon their construction, the evaluation time of the covariance matrices is approximately linear w.r.t. the number of test points. On the other hand, Laplace-Beltrami based method require the one-time eigenfunction decomposition—which is time demanding—but then provides an overall faster inference.

Our method shows its maximum strength in all cases in which the mesh is not known upfront, e.g., for time-varying surfaces, it is small to medium sized (up to 5\,000 faces), and it is limited to the computation of the shortest path on \mathcal{M}, which is delegated to external libraries.

V-C Mesh Flow Matching

To fairly compare our proposed Mesh Flow Matching against the baseline, we integrated our mesh operators in the official implementation of the RFM framework provided by [4]. As in the reference paper, the surface used for the experiment is the Stanford Bunny, decimated to 5\,000 faces, and Spot the Cow. To construct the non-trivial target distribution that must be learned from the generative model, 500\,000 samples are drawn from the mesh by computing the i-th eigenfunction of the Laplace-Beltrami operator (26) and sampling points on \mathcal{M} proportionally to the thresholded eigenfunction value. The dataset is split into training, validation, and test sets using an 80–10–10 ratio. An example of target data distribution and sample reconstruction with our mesh RFM method is reported in Fig. 12.

To evaluate the regression performance, we compute the negative log-likelihood (NLL) by averaging the log density \log\pi_{\textrm{data}}(p) of all samples in the test dataset. Following [34], log-density is computed as

\log\pi_{\textrm{data}}(p)=\log\pi_{\textrm{ref}}(p_{0})-f_{0}(p),

where f_{t} is obtained solving backward in time from t=1 to t=0 the ODE

\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}p_{t}\\
f_{t}(p)\end{pmatrix}=\begin{pmatrix}v_{t}(p_{t})\\
-\textrm{div}_{\mathfrak{g}}(v_{t})(p_{t})\end{pmatrix}

subject to the initial condition \big(p_{1},f_{1}(p)\big)^{\top}=(p,0)^{\top}; here, \textrm{div}_{\mathfrak{g}} is the Riemannian divergence.

Table IV reports the test NLL—evaluated at the model checkpoint with lowest validation NLL—of our method, and the RFM baseline which uses diffusion and biharmonic premetrics computed on the first 200 eigenfunctions. All models have been trained with the same parameters and seeds, using a Dense NN with 6 layers of 512 neurons each with the Swish activation function [55], and using the ADAM optimiser with a learning rate of 10^{-4}. Across all experiments, our method consistently achieves lower test NLL, indicating improved generative performance when using geometry-aware operators defined directly on the mesh.

Computation performance

The claimed advantage of Laplace-Beltrami based spectral operators, is that they limit the computation time by only considering a fixed, relatively small number of eigenfunctions—200 in this case. However, since the integration of the spectral measures doesn’t guarantee the path to lie on the mesh, closest-point projections are required by the baseline algorithm at every step, effectively slowing down the integration process. On the other hand, our method—by using the exponential map (14)—is always guaranteed, during ODE integration, to generate paths which are consistent with the underlying surface. Additionally, during the training (Algorithm 4), the time-intensive geodesic path computation between samples p_{0} and p_{1} is amortised by the cheap upsampling over time of the vector field u_{t}.

From our experience, performing 100\,000 update steps of the network with batches of 128 samples, and calling a validation routine every 2\,000 steps to save model checkpoints, typically takes 1 day and 6 hours for the baseline method, while the proposed method takes on average 20 hours.

VI Experiments

VI-A Handwritten symbols learning and generalisation

The key feature of DMPs is their high generalisation capability, which enables the execution of motion policies to new, unseen initial and final condition while preserving the topology of the learned path by exploiting spatial and temporal invariance properties. As highlighted in (22) and (23), exploiting the definition of surface operators on discrete manifolds, it is possible to learn motion policies that can be transferred to different meshes. We show this with an experiment in which user-acquired demonstration are generalised to different surfaces to enable a robotic manipulator drawing on them.

For this experiment, we created 4 surfaces generated as the plot of continuous functions of type H:\mathbb{R}^{2}\rightarrow\mathbb{R}; upon fixing a triangulated grid on the 20\times 20\,cm squared domain, we construct the discrete approximation of the surface (x,y)\mapsto\big(x,y,H(x,y)\big), and we 3D printed these samples. To recreate a realistic scenario of a smooth surface, the 3D printed components is generated with a high polygonal density, with the triangulated grid having 500\times 500 points; during experiments, instead, we relied on a coarser mesh defined on grids ranging from 40\times 40 to 80\times 80 points. To enable robot motion on the surface, the mesh samples have been localised relative to the base of the Franka Research 3 manipulator [56] that we used as robotic platform. To compensate localisation inaccuracies, we employed an hybrid force-impedance control; particularly, the robot was prompted to track a 5\,N force along the pencil direction, while other axis were set on impedance control with linear stiffness of 4\,000\,N/m, rotational stiffness 100\,Nm/rad, and critical damping.

To ease robot teaching, we designed a simple and intuitive web-based Graphical User Interface (GUI) in which non-expert users can sketch shapes on a tablet, and then decide which type of DMP to use through buttons. The system immediately learns the flat demonstration and generalises it to all target meshes, creating 40\,000 trajectory samples typically in less then 5\,s; the GUI also enables updating DMPs parameters—such as starting and final configuration, or motion centre—and control the robot behaviour, like starting or stopping the execution of the motion on a specific surface. Since the proposed DMPs on surfaces provide only position information for the robots end effector, we select the orientation of the pencil to be normal to the surface by exploiting the normal map given by the mesh itself. This choice leads to non-differentiable orientation trajectories when transitioning between faces; for this reason, we firstly precompute the whole trajectory, and then use a symmetric sliding-window filter to smoothen the orientation component.

Figure 13 effectively shows how the proposed operators on meshes, paired with DMPs (22) and (23), can generalise the motion from the flat plane to the different types of surfaces. Since the operators require no smoothness condition of the mesh \mathcal{M}, our method works also when the underlying function H is non-differentiable in some point locus, as in the case depicted in the last column of that figure. We refer to the accompanying video for more experiments.

VI-B Surface polishing

In this experiment, we show how our methodology can be easily integrated in a larger framework for fully autonomous surface polishing, where the surface is reconstructed from sensory data. We consider the setup in Fig. 14, where the robot is equipped with an Intel Realsense D405 RGB-D camera. The depth camera is used to generate a point cloud of the observed scene, which enables online surface reconstruction and, through simple computer vision algorithms, the generation of polishing motions on the target surface using our rhythmic DMP.

As target workpiece, we use a blue-coloured component with non-flat surface, and draw a line on it with a orange chalk; by using range thresholding on the HSL-channels of the colour image, we masked the drawn line and proceeded fitting a polynomial on it, as shown in Fig. 14(a). To improve mesh reconstruction reliability, we use a waterfront expanding algorithm to segment the full point cloud in the region of interested given by the polynomial path. To further reduce computational burden, we used Voxel downsampling with cell size of 5\,mm on the segmented point cloud, and then proceeded in using the Poisson surface reconstruction algorithm [57] to generate a watertight, manifold mesh that we can use for planning. In our experiments, the proposed pipeline takes from 1 to 10 seconds, based on point-cloud quality; for the experiment shown in Fig. 14, the generated mesh had 1\,344 vertices and 2\,617 faces.

To achieve the polishing behaviour over the surface, we first trained a rhythmic DMP (23) on a circular motion on the flat surface. We then generalised the policy to the point cloud-reconstructed mesh, and achieve a polishing-like behaviour by slowly shifting the DMP centre g along the properly projected fitted polynomial. Fixed the velocity of the shifting centre along the path to 5\,mm/s and a time-per-revolution of 2\,s, the algorithm generates offline the full position trajectory at 100\,Hz typically in less then five seconds. The resulting motion is overlaid on the acquired point cloud in Fig. 14(b).

For robotic execution, we equipped the end-effector with a soft sponge, and leverage the hybrid force-impedance controller to exert a constant force on the surface. The generated trajectory is upsampled to 1\,kHz and, as in Sec. VI-A, we assign the orientation to be normal to the surface; to ensure smoothness of the orientation signal, we used a symmetric sliding window of 1\,s. Figure 14(c)–(e) reports the surface before and after the robot polishing, showing how the proposed methodology can be effectively integrated in real tasks where the object surface is not known upfront.

VII Discussion and Conclusion

In this work, we developed a suite of geometrically consistent operators—logarithmic and exponential map, parallel transport, and ambient space projections—that extend Riemannian differential geometry to the realm of discrete surfaces modelled by polyhedral meshes. The core of our approach is to work directly on the intrinsic geometry of the surfaces, rather than considering the surfaces as mere constraints to be enforced a posteriori. We successfully applied these operators to a diverse set of geometry-aware learning methods relevant to the robotics community, namely DMPs, GPs, and RFM.

Concerning DMPs, by encoding the forcing term within a fixed tangent cone and reconstructing it through parallel transport, we eliminate the limit-cycle instability of prior approaches, constructing motion policies that generalise better across multiple surfaces. We validated our approach in simulation—comparing against the MeshDMP framework [9]—and on a real end-to-end LfD pipeline, with a Franka Research 3 manipulator reproducing tablet-drawn trajectories. Both the simulation and the real-world experiment demonstrate the numerical efficiency of the approach. As regards GPs, replacing spectral Laplace-Beltrami measures [10] with a geodesic-based kernel removes all smoothness assumptions for the underlying surface, enabling regression at arbitrary locations on the surface without interpolation. This technique is particularly effective on coarse meshes. Finally, for RFMs, the integration of mesh-specific operators improves generative quality compared to spectral baselines, whilst reducing training time by approximately half.

Taken together, these results suggest a remarkable advance towards bridging geometry-aware learning theory and the mesh representations employed in robotics. The cornerstone of our construction is a principled application of discrete differential geometry.

Limitations

Although we believe that the advances brought about by our approach are substantial, some of the present limitations should not be overlooked, as they also offer interesting opportunities for future investigations. First, the geodesic pre-computation remains a bottleneck. Since it relies on algorithms with complexity \mathcal{O}(N^{2}\log N), where N is the number of faces of the mesh, this computation becomes a significant barrier when the surface exceeds 10\,000 faces.

Second, the geodesic-based kernel cannot be guaranteed to be positive-definite when the lengthscale exceeds some critical value \kappa_{th}. This unavoidable problem comes with the risk of generating inconsistent GPs when the kernel is evaluated beyond the critical value. While Algorithm 3 offers a practical safeguard, a better characterisation of the threshold \kappa_{th} is required—particularly by relating this coefficient to surface curvature and topology, as well as the number of sampled points.

Future directions

Whilst this paper presents a strong mathematical foundation and three important application paradigms, we believe that much remains to be done to extend the framework to a broader set of applications.

One promising research avenue is to tackle the challenges posed by the extension of the formulation to learn force trajectories on surfaces and, at the same time, to transfer them across surfaces with different topologies. Another open research direction is to extend the framework to time-varying surfaces, e.g., to enable motion planning for surgical applications on human bodies whose pose is estimated online via parameterised meshes that evolve over time with the patient’s movements.

References

1. [1] L. Sentis and O. Khatib (2005) Synthesis of whole-body behaviors through hierarchical control of behavioral primitives. International Journal of Humanoid Robotics 2 (4), pp. 505–518 (en).
2. [2] S. Calinon (2020) Gaussians on Riemannian manifolds: applications for robot learning and adaptive control. IEEE Robotics and Automation Magazine (RAM) 27 (2), pp. 33–45. External Links: Document
3. [3] N. Jaquier, L. Rozo, D. G. Caldwell, and S. Calinon (2021) Geometry-aware manipulability learning, tracking, and transfer. International Journal of Robotics Research (IJRR) 40 (2-3), pp. 624–650 (en).
4. [4] R. T. Q. Chen and Y. Lipman (2024) Flow matching on general geometries. In 12th International Conference on Learning Representations (ICLR), External Links: 2302.03660
5. [5] F. J. Abu-Dakka and V. Kyrki (2020) Geometry-aware dynamic movement primitives. In 2020 IEEE International Conference on Robotics and Automation (ICRA), Vol. , pp. 4421–4426. External Links: Document
6. [6] M. Berger, A. Tagliasacchi, L. M. Seversky, P. Alliez, G. Guennebaud, J. A. Levine, A. Sharf, and C. T. Silva (2017) A survey of surface reconstruction from point clouds. Computer Graphics Forum 36 (1), pp. 301–329.
7. [7] A. Hubeli and M. Gross (2000) A survey of surface representations for geometric modeling. Technical report ETH Zurich, Department of Computer Science, ETH Zurich, Department of Computer Science.
8. [8] M. Pantic, L. Ott, C. Cadena, R. Siegwart, and J. Nieto (2021) Mesh manifold based riemannian motion planning for omnidirectional micro aerial vehicles. IEEE Robotics and Automation Letters 6 (3), pp. 4790–4797. External Links: Document
9. [9] M. Dalle Vedove, F. J. Abu-Dakka, L. Palopoli, D. Fontanelli, and M. Saveriano (2025) MeshDMP: motion planning on discrete manifolds using dynamic movement primitives. In 2025 IEEE International Conference on Robotics and Automation (ICRA), Vol. , pp. 895–901. External Links: Document
10. [10] V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth (2020) Matérn gaussian processes on riemannian manifolds. In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS ’20, Red Hook, NY, USA. External Links: ISBN 9781713829546
11. [11] A. I. Bobenko and Y. B. Suris (2023) Discrete differential geometry: integrable structure. American Mathematical Society, Providence, RI (en).
12. [12] K. Polthier and M. Schmies (2006) Straightest geodesics on polyhedral surfaces. In ACM SIGGRAPH 2006 Courses on - SIGGRAPH ’06, SIGGRAPH ’06, pp. 30. External Links: Document
13. [13] X. Xue, J. Dong, Z. Lu, and N. Wang (2023) A robotic learning and generalization framework for curved surface based on modified dmp. Robotics and Autonomous Systems 160, pp. 104323.
14. [14] A. G. Billard, S. Calinon, and R. Dillmann (2016) Learning from humans. Springer handbook of robotics, pp. 1995–2014.
15. [15] M. Saveriano, F. J. Abu-Dakka, A. Kramberger, and L. Peternel (2023) Dynamic movement primitives in robotics: a tutorial survey. The International Journal of Robotics Research 42 (13), pp. 1133–1184. External Links: ISSN 0278-3649
16. [16] A. Paraschos, C. Daniel, J. R. Peters, and G. Neumann (2013) Probabilistic movement primitives. Advances in neural information processing systems 26.
17. [17] Y. Huang, L. Rozo, J. Silvério, and D. G. Caldwell (2019) Kernelized movement primitives. The International Journal of Robotics Research 38 (7), pp. 833–852.
18. [18] S. Calinon (2016) A tutorial on task-parameterized movement learning and retrieval. Intelligent service robotics 9, pp. 1–29.
19. [19] A. Ude, B. Nemec, T. Petrić, and J. Morimoto (2014) Orientation in cartesian space dynamic movement primitives. In 2014 IEEE International Conference on Robotics and Automation (ICRA), Vol. , pp. 2997–3004. External Links: Document
20. [20] F. J. Abu-Dakka and M. Saveriano (2020) Variable impedance control and learning-a review. Front. Robot. AI 7, pp. 590681 (en).
21. [21] L. Koutras and Z. Doulgeri (2020) A correct formulation for the orientation dynamic movement primitives for robot control in the cartesian space. In Proceedings of the Conference on Robot Learning, L. P. Kaelbling, D. Kragic, and K. Sugiura (Eds.), Proceedings of Machine Learning Research, Vol. 100, pp. 293–302. External Links: Link
22. [22] J. Kim, M. Sung, Y. Choi, J. Park, and W. K. Chung (2025) Impedance control design framework using commutative map between SE(3) and \mathfrak{se}(3). IEEE Transactions on Robotics 41 (), pp. 6193–6212. External Links: Document
23. [23] F. J. Abu-Dakka, M. Saveriano, and V. Kyrki (2024) A unified formulation of geometry-aware discrete dynamic movement primitives. Neurocomputing 598, pp. 128056. External Links: ISSN 0925-2312, Document
24. [24] N. D. Ratliff, J. Issac, D. Kappler, S. Birchfield, and D. Fox (2018) Riemannian motion policies. External Links: arXiv:1801.02854
25. [25] E. Dimas and D. Briassoulis (1999) 3D geometric modelling based on NURBS: a review. Adv. Eng. Softw. 30 (9-11), pp. 741–751 (en).
26. [26] I. Boier-Martin and D. Zorin (2004) Differentiable parameterization of catmull-clark subdivision surfaces. In ACM Symposium on Geometry Processing, pp. 155–164. External Links: ISBN 3905673134, Link, Document
27. [27] G. P. T. Choi, Y. Leung-Liu, X. Gu, and L. M. Lui (2020) Parallelizable global conformal parameterization of simply-connected surfaces via partial welding. SIAM J. Imaging Sci. 13 (3), pp. 1049–1083 (en).
28. [28] B. Vallet and B. Lévy (2008) Spectral geometry processing with manifold harmonics. Comput. Graph. Forum 27 (2), pp. 251–260 (en).
29. [29] K. Crane, F. de Goes, M. Desbrun, and P. Schröder (2013) Digital geometry processing with discrete exterior calculus. In ACM SIGGRAPH 2013 CoursesSIGGRAPH ’13: Special Interest Group on Computer Graphics and Interactive Techniques Conference, New York, NY, USA.
30. [30] G. Patané (2014) Laplacian spectral distances and kernels on 3d shapes. Pattern Recognition Letters 47, pp. 102–110. Note: Advances in Mathematical Morphology External Links: ISSN 0167-8655, Document, Link
31. [31] L. Koestler, D. Grittner, M. Moeller, D. Cremers, and Z. Lähner (2022) Intrinsic neural fields: learning functions on manifolds. In European Conference on Computer Vision, S. Avidan, G. Brostow, M. Cissé, G. M. Farinella, and T. Hassner (Eds.), Cham, pp. 622–639.
32. [32] A. Lang and M. Pereira (2023) Galerkin–Chebyshev approximation of gaussian random fields on compact riemannian manifolds. BIT 63 (4) (en).
33. [33] A. Feragen, F. Lauze, and S. Hauberg (2015) Geodesic exponential kernels: when curvature and linearity conflict. In IEEE Conference on Computer Vision and Pattern Recognition, pp. 3032–3042. External Links: Document
34. [34] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud (2018) Neural ordinary differential equations. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, Red Hook, NY, USA, pp. 6572–6583.
35. [35] Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le (2023) Flow matching for generative modeling. In The Eleventh International Conference on Learning Representations,
36. [36] X. Liu, C. Gong, and Q. Liu (2023) Flow straight and fast: learning to generate and transfer data with rectified flow. In The Eleventh International Conference on Learning Representations,
37. [37] M. P. do Carmo (2013) Riemannian geometry. 14., corrected print. edition, Mathematics: theory & applications, Birkhäuser, Boston. Note: Aus dem Portug. übers. - Literaturverz. S. [292] - 296 External Links: ISBN 3764334908
38. [38] W. Hassler (1944) The self-intersections of a smooth n-manifold in 2n-space. Ann. Math. 45 (2), pp. 220.
39. [39] L. De Floriani and A. Hui (2005) Data structures for simplicial complexes: an analysis and a comparison. In Proceedings of the Third Eurographics Symposium on Geometry Processing, SGP ’05, Goslar, DEU, pp. 119–es. External Links: ISBN 390567324X
40. [40] A. J. Ijspeert, J. Nakanishi, H. Hoffmann, P. Pastor, and S. Schaal (2013) Dynamical movement primitives: learning attractor models for motor behaviors. Neural Computation 25 (2), pp. 328–373. External Links: Document
41. [41] C. E. Rasmussen and C. K. I. Williams (2005) Gaussian processes for machine learning. The MIT Press. External Links: ISBN 9780262256834, Document
42. [42] F. Clarke (2013) Functional analysis, calculus of variations and optimal control. Springer London. External Links: ISBN 9781447148203, Document, ISSN 2197-5612
43. [43] J. S. B. Mitchell, D. M. Mount, and C. H. Papadimitriou (1987) The discrete geodesic problem. SIAM Journal on Computing 16 (4), pp. 647–668. External Links: ISSN 1095-7111, Document
44. [44] J. Chen and Y. Han (1990) Shortest paths on a polyhedron. In Proceedings of the sixth annual symposium on Computational geometry - SCG ’90, SCG ’90, pp. 360–369. External Links: Document
45. [45] S. Xin and G. Wang (2009) Improving chen and han’s algorithm on the discrete geodesic problem. ACM Transactions on Graphics 28 (4), pp. 1–8. External Links: ISSN 1557-7368, Document
46. [46] M. P. do Carmo (2016) Differential geometry of curves and surfaces. Revised & updated second edition edition, Mathematics, Dover Publications,Inc., Mineola, New York. External Links: ISBN 9780486817972
47. [47] M. Kniely and W. Ring (2014) Riemannian methods for optimization in a shape space of triangular meshes. Inverse Problems in Science and Engineering 23 (6), pp. 1011–1039. External Links: ISSN 1741-5985, Document
48. [48] C. Ericson (2004) Real-Time collision detection. Morgan Kaufmann series in interactive 3D technology, CRC Press, London, England (en).
49. [49] The CGAL Project (2026) CGAL user and reference manual. 6.1.1 edition, CGAL Editorial Board. External Links: Link
50. [50] S. Kiazyk, S. Loriot, and É. C. de Verdière (2026) Triangulated surface mesh shortest paths. In CGAL User and Reference Manual, External Links: Link
51. [51] K. Crane, U. Pinkall, and P. Schröder (2013) Robust fairing via conformal curvature flow. ACM Transactions on Graphics (TOG) 32 (4), pp. 1–10.
52. [52] B. Curless and M. Levoy (1996) A volumetric method for building complex models from range images. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniquesInternational Conference on Computer Graphics and Interactive Techniques,
53. [53] P. Whittle (1963) Stochastic processes in serveral dimensions. Bulletin of the International Statistical Institute 40, pp. 974–994.
54. [54] P. Mostowsky, V. Dutordoir, I. Azangulov, N. Jaquier, M. J. Hutchinson, A. Ravuri, L. Rozo, A. Terenin, and V. Borovitskiy (2024) The geometrickernels package: heat and matérn kernels for geometric learning on manifolds, meshes, and graphs. arXiv:2407.08086.
55. [55] P. Ramachandran, B. Zoph, and Q. V. Le (2018) Searching for activation functions. In International Conference on Learning Representations,
56. [56] S. Haddadin, S. Parusel, L. Johannsmeier, S. Golz, S. Gabl, F. Walch, M. Sabaghian, C. Jähne, L. Hausperger, and S. Haddadin (2022) The franka emika robot: a reference platform for robotics research and education. IEEE Robotics & Automation Magazine 29 (2), pp. 46–64.
57. [57] M. Kazhdan, M. Bolitho, and H. Hoppe (2006) Poisson surface reconstruction. In Proceedings of the Fourth Eurographics Symposium on Geometry Processing, SGP ’06, Goslar, DEU, pp. 61–70. External Links: ISBN 3905673363