Tactile Ergodic Coverage on Curved Surfaces

1 Tactile Ergodic Coverage on Curved Surfaces Cem Bilaloglu∗, Tobias Löw∗, and Sylvain Calinon Abstract—In this article, we present a feedback control method for tactile coverage tasks such as cleaning or surface inspection. Although these tasks are challenging to plan due to the complexity of continuous physical interactions, the coverage target and progress can be effectively measured using a camera and encoded in a point cloud. We propose an ergodic coverage method that operates directly on point clouds, guiding the robot to spend more time on regions requiring more coverage. For robot control and contact behavior, we use geometric algebra to formulate a task-space impedance controller that tracks a line while simultaneously exerting a desired force along that line. We evaluate the performance of our method in kinematic simulations and demonstrate its applicability in real-world experiments on kitchenware. Our source codes, experiment videos, and data are available at https://sites.google.com/view/tactile-ergodic-control/. Index Terms—Tactile Robotics, Ergodic Coverage, Geometric Algebra I. INTRODUCTION The long-term vision of robotics is to assist humans with daily tasks. The success of robot vacuum cleaners and lawnmowers as consumer products highlights the potential of robotic assistance for common household chores [1]. These tasks involve covering a region in a repetitive and exhaustive manner. Currently, these robots are limited to relatively large, planar surfaces, and even navigating slopes remains challenging [2], [3]. Other daily tasks, such as washing dishes or grocery items, present even greater challenges due to the complex physical interactions with intricate, curved surfaces. Similarly, numerous coverage tasks on curved surfaces arise in industrial and medical applications. In industrial settings, such tasks include surface operations that remove material, such as sanding [4], polishing [5], [6] or deburring [7] as well as surface inspection tasks leveraging contact [8]. In medical settings, similar applications range from mechanical palpation [9], [10] and ultrasound imaging [11], [12] to massage [13], [14] and bed bathing [15], [16]. Last but not least, datasets combining the tactile properties of objects with their shape and visual appearance remain scarce and expensive to collect, as they rely on teleoperation [17]. Thus, tactile coverage is critical for automating the collection of tactile datasets that complement visual ones. The problem definitions of this diverse range of settings and applications can be distilled into two key requirements: (i) tactile interactions with This work was supported by the State Secretariat for Education, Research and Innovation in Switzerland for participation in the European Commission’s Horizon Europe Program through the INTELLIMAN project (https: //intelliman-project.eu/, HORIZON-CL4-Digital-Emerging Grant 101070136) and the SESTOSENSO project (http://sestosenso.eu/, HORIZON-CL4Digital-Emerging Grant 101070310). ∗Equal contribution. The authors are with the Idiap Research Institute, Martigny, Switzerland and with the Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland. [email protected]; [email protected]; [email protected] Diffusion Coverage Gradient field Measure coverage Update target Coverage Loop Vision Tactile Control Force Fig. 1: Overview of our feedback control method for tactile coverage. Left: We measure the surface and the red target using the camera and encode them in a point cloud. Bottom-right: We diffuse the target and use its gradient field to guide the coverage. Then, we close the loop by measuring the actual coverage with the camera and use it as the next target. Topright: We measure the tactile interaction forces using the force sensor and the tool orientation using the joint positions. We solve the geometric task-space impedance control problem using a line target and a force target along the line. a possibly non-planar surface and (ii) a continuous trajectory of contact points covering a region of interest on the surface. Accordingly, this article addresses the overarching problem of tactile coverage on curved surfaces. Tactile tasks involve multiple contact interactions with the environment, making these systems notoriously difficult to control [18]. While humans solve these tasks effortlessly, they remain extremely challenging for robots. For instance, when cleaning an object, achieving adequate coverage depends on recognizing dirt, understanding the object’s material, and assessing their interaction to determine the required contact force for removal. Consequently, the success of coverage depends on unknown or difficult-to-measure parameters, making it challenging to model all interactions. Without an accurate model, motion planning is prone to failure. By analyzing previous research [19] and observing how humans address these challenges, we argue that humans bypass the complexities of planning by solving the simpler closed-loop control problem. Humans leverage visual and tactile feedback for online adaptation. Similarly, robots can measure progress in tactile coverage tasks using vision, turning the task of identifying uncovered regions into an image segmentation problem which has been addressed using various model-based [20], [21] or learningbased algorithms [16], [22]. However, determining how to control a robot to cover these target regions on curved surfaces remains an open challenge. arXiv:2402.04862v3 [cs.RO] 31 Mar 2025 2 Existing research on coverage has primarily focused on coverage path planning, which involves optimizing a path to ensure that a specified region of interest is covered within a set time frame. Traditionally, the underlying assumption is that visiting each point in the region of interest only once is sufficient for full coverage, an assumption that is reasonable for simple interactions, but not for many tactile tasks. Tactile interactions are often too complex to model deterministically, making it challenging to ensure full coverage after a single visit. Instead, for a cleaning task, a relatively dirty region requires more visits compared to a less dirty region. Similarly, in a surface inspection task, regions requiring higher precision demand more visits to compensate for sensor uncertainty. Furthermore, the robot is expected to keep in contact with the surface while moving, which significantly increases the cost of movement. This cost depends on the geodesic distance on the surface rather than the Euclidean distance. Therefore, naive sampling strategies that fail to account for the cost or constraints of movement and/or surface geometry are unsuitable for tactile coverage tasks. In contrast, ergodic coverage [23] controls the trajectories of dynamical systems by correlating the average time spent in a region to the target spatial distribution. Therefore, ergodic coverage incorporates the motion model as the system dynamics and directly controls the coverage trajectories by using the spatial distribution measured by the vision system. Considering these challenges, we present a closed-loop tactile ergodic control method that operates on point clouds for tactile coverage tasks. Using point clouds enables us to acquire the target object and spatial distribution at runtime using vision, measure coverage progress, and compensate for unmodeled dynamics in tactile coverage tasks. Our method then constrains the ergodic control problem to arbitrary surfaces to cover a target spatial distribution on the surface. We propagate coverage information by solving the diffusion equation on point clouds, which we compute in real-time by exploiting the surface’s intrinsic basis functions called Laplacian eigenfunctions. These eigenfunctions generalize the Fourier series to manifolds (i.e., curved spaces). In order to exert a desired force on the surface while moving, we formulate a geometric task-space impedance controller using geometric algebra. This controller uses surface information to track a line target that is orthogonal to the surface while simultaneously exerting the desired force in the direction of that line. Notably, the geometric formulation ensures that these two objectives do not conflict with each other and can therefore be included in the same control loop without requiring exhaustive parameter tuning. In summary, our proposed closed-loop tactile ergodic control method offers the following contributions: •formulating the tactile coverage as closed-loop ergodic control problem on curved surfaces; •closing the coverage loop by solving ergodic control problem on point clouds using diffusion; •achieving real-time frequencies by computing the diffusion using Laplacian eigenfunctions; •contact line and force tracking without conflicting objectives. The rest of the article is organized as follows. Section II describes related work. Section III provides the mathematical background. Section IV presents our method. In Section V, we demonstrate the effectiveness of our method in simulated and real-world experiments. Finally, we discuss our results in Section VI. II. RELATED WORK The majority of the coverage methods consider the problem from a planning perspective and are generally known as coverage path planning (CPP) algorithms [24]–[26]. Although these methods can handle planar regions with various boundaries [27]–[29], their extension to curved surfaces imposes limiting assumptions, such as projectively planar [30] or pseudoextruded surfaces [31]. Additionally, CPP methods assume that the coverage target is uniformly distributed in space. Extending CPP methods to account for spatial correlations in the information leads to informative path planning (IPP) [32]. Most IPP and CPP approaches address a variant of the NP-hard traveling salesman problem [33], which limits their scalability as domain complexity increases. Consequently, existing methods are either open-loop [34] or impose limiting assumptions, such as convexity, for online planning updates [32]. Closely related to coverage is the problem of exploration, where the environment is initially unknown, and robots gather information using onboard sensors [35], [36]. Tactile exploration is particularly necessary for gathering information on surfaces that can only be acquired through contact [37]. A notable example is non-invasive probing (palpation) of tissue stiffness, which aids in disease diagnosis or surgery by providing additional anatomical information. For this purpose, Gaussian processes (GP) have been used for discrete [9] and continuous [38] probing to map tissue stiffness. While GPbased approaches effectively guide sampling locations, they do not account for the robot’s dynamics. This limitation was later addressed by using trajectory optimization to actively search for tissue abnormalities [39]. Unlike other sensing modalities that depend solely on position, tactile interactions also depend on conditions such as relative velocity and contact pressure [40]. To address this, methods have been developed to model forces [41] and more complex interactions between robotic tools and surfaces [42]. The complexity of the problem increases further if we consider scenarios with a robot physically interacting with the environment. For example, in tasks like surface finishing (e.g., polishing, sanding, grinding), the surface itself changes, as material is removed [43]. Similarly, in cleaning tasks, the robot’s actions affect the distribution of dirt on the surface [44]. To avoid complex modeling, there are approaches either relying on reinforcement learning [45] or deep learning [46]. In a very similar setting to ours, a manipulator was used to clean the stains on a curved surface by performing multiple passes [20]. However, this work used a sampling-based planner, which required to predefine the maximum number of cleaning passes. In contrast, we relate the target distribution (e.g., stain) directly to feedback control without requiring any task-specific assumptions. 3 In tactile coverage scenarios, visiting a region once can not guarantee full coverage, and predicting how many times the robot should revisit a particular spot is challenging. Consequently, defining a time horizon for trajectory optimization is difficult, as the quality of the result would be significantly affected by this hard-to-make choice. Instead, ergodic control relates how often the robot should revisit a particular spot to the target density at that spot. In this context, ergodic describes a dynamical system in which the time averages of functions along its trajectories are equal to their spatial averages [47]. The key advantage of ergodic control is its ability to handle arbitrary spatial target distributions without requiring a predefined time horizon. When the spatial target distribution is measurable, the ergodic controller can use this feedback to direct the system to visit regions with higher spatial probabilities more frequently. Recent findings have demonstrated that ergodicity is not merely a heuristic [48]; it is the optimal method for collecting independent and identically distributed data while accounting for system dynamics. Ergodic control was introduced in the seminal work by Mathew and Mezi´ c [23], which presented the spectral multiscale coverage (SMC) algorithm. SMC is a feedback control law based on the Fourier decomposition of the target distribution and robot trajectories, where multiscale aspect prioritizes low-frequency components over high-frequency ones, corresponding to starting with large-scale spatial motions before refining finer details. Since this behavior is achieved through a myopic feedback controller rather than an offline planner, the ergodic controller remains effective even when motion is obstructed [49]. Recent works have adapted SMC’s objective within a trajectory optimization framework to incorporate additional objectives, such as obstacle avoidance [50], timeoptimality [51] and energy-awareness [52]. Ergodic control has been used for tactile coverage and exploration in applications such as non-parametric shape estimation [53] and table cleaning through learning from demonstration [54]. However, all these formulations, which rely on the Fourier decompositionbased ergodic metric, are limited to rectangular domains in Euclidean space. The first attempt to extend the ergodic control to Riemannian manifolds [55] utilized Laplacian eigenfunctions, which generalize the Fourier series to curved spaces. However, this approach was restricted to homogeneous manifolds, such as spheres and tori, where closed-form expressions for the Laplacian eigenfunctions are available. More recently, the kernel ergodic metric [56] was introduced as an alternative to SMC’s ergodic metric, enabling extensions to Lie groups and offering improved computational scalability. Nonetheless, arbitrary curved surfaces collected using sensors, such as point clouds, lack both the group structure and the homogeneous manifold properties, presenting additional challenges. Another alternative to SMC is the heat equation-driven area coverage (HEDAC) algorithm [57], which uses the diffusion equation, a second-order partial differential equation (PDE), to propagate information about uncovered regions to agents across the domain. Similar to SMC, the original HEDAC implementation was restricted to rectangular domains and lacked collision avoidance. Subsequent extensions have adapted HEDAC to planar meshes with obstacles [58], maze exploration [59], and CPP on non-planar meshes [60]. However, its application on curved surfaces remains limited to meshes and offline planning due to the heavy pre-processing required. In addition to its use in HEDAC, the diffusion equation is widely used in geometry processing tasks, ranging from geodesic computation [61] to learning on surfaces [62]. Its key advantage lies in its ability to account for surface geometry while remaining agnostic to the underlying representation and discretization [62]. The diffusion equation is governed by a second-order differential operator called the Laplacian which can be computed for arbitrary surfaces represented as meshes or point clouds using various discretization schemes [63]–[65]. In this work, we use a recent approach proposed by Sharp et al. which provides a robust and efficient implementation [66], capable of handling partial and noisy point clouds. III. BACKGROUND A. Ergodic Control using Diffusion The ergodic control objective correlates the time that a coverage agent spends in a region to the probability density specified in that region. The HEDAC method [57] encodes the coverage objective in the domain x∈Ωat time tusing a virtual source term s(x, t) = max p(x)−c(x, t),02,(1) where p(x)is the probability distribution corresponding to the coverage target and c(x, t)is the normalized coverage of the Nvirtual coverage agents over the domain c(x, t) = ˜c(x, t) RΩ˜c(x, t)dx.(2) A single agent’s coverage is the convolution of its footprint φ(r)with its trajectory xi(t′). Then, the total coverage becomes the time-averaged sum of these convolutions ˜c(x, t) = 1 Nt N X i=1 Zt 0 φx−xi(t′)dt′.(3) HEDAC diffuses the source term across the domain Ωand computes the potential field u(x, t)using the stationary (˙u(x, t) = 0) diffusion α∆u(x, t)−u(x, t) + s(x, t)=0,(4) with the diffusion coefficient α > 0and the Laplacian operator ∆. The Laplacian is a second-order differential operator which reduces to the sum of the second partial derivatives in Euclidean spaces ∆f=∇·∇f= n X i=1 ∂2f ∂x2 i ,∀x∈Rn.(5) The stationary diffusion (4) governs the potential field within the interior of the domain Ω, while the behavior on the boundary ∂Ωis dictated by the zero-Neumann boundary condition n· ∇u(x, t) = 0,∀x∈∂Ω,(6) 4 where nrepresents the outward unit normal vector to the boundary ∂Ω. To guide the i-th coverage agent, HEDAC utilizes the smooth gradient field of the diffused potential u(x, t)and simulates first-order dynamics [67] ˙ xi=∇u(xi, t).(7) B. Conformal Geometric Algebra Here, we introduce conformal geometric algebra (CGA) with a focus on the mathematical background necessary to understand the methods used in this article. We will use the following notation throughout the paper: xto denote scalars, xfor vectors, Xfor matrices, Xfor multivectors and Xfor matrices of multivectors. The inherent algebraic product of geometric algebra is called the geometric product ab =a·b+a∧b,(8) which (for vectors) is the sum of an inner ·and an outer ∧ product. The inner product is the metric product and therefore depends on the metric of the underlying vector space over which the geometric algebra is built. The underlying vector space of CGA is R4,1, which means there are four basis vectors squaring to 1 and one to -1. The outer product, on the other hand, is a spanning operation that effectively makes subspaces of the vector space elements of computation. These subspaces are called blades. In the case of CGA, there are 32 basis blades of grades 0 to 5. The term grade refers to the number of basis vectors in a blade that are factorizable under the outer product. Vectors, consequently, are of grade 1 and the outer product of two independent vectors, called bivectors, are of grade 2. A general element of geometric algebra is called a multivector. In practice, CGA actually applies a change of basis by introducing the two null vectors e0and e∞, which can be thought of as a point at the origin and at infinity, respectively. Since the Euclidean space is embedded in CGA, we can embed Euclidean points xto conformal points Pvia the conformal embedding P=C(x) = e0+x+1 2x2e∞.(9) In general, geometric primitives in geometric algebra are defined as nullspaces of either the inner or the outer product, which are dual to each other. The outer product nullspace (OPNS) is defined as NOG(X) = x∈R3:C(x)∧X= 0.(10) A similar expression can be found for the inner product nullspace. The conformal points are the basic building blocks to construct other geometric primitives in their OPNS representation. The relevant primitives for this work are lines L=P1∧P2∧e∞,(11) which can be constructed from two points and a point at infinity, planes E=P1∧P2∧P3∧e∞,(12) which can be constructed from three points and a point at infinity and spheres S=P1∧P2∧P3∧P4,(13) which can be constructed from four points. Rigid body transformations in CGA are achieved using motors M, which are exponential mappings of dual lines, i.e. bivectors (essentially, the screw axis of the motion). Note that motors can be used to transform any object in the algebra, i.e. they can directly be used to transform the previously introduced points, lines, planes and spheres, by a sandwiching operation X′=MX f M, (14) where is f Mis the reverse of a motor. The forward kinematics of serial kinematic chains can be found as the product of motors, i.e. M(q) = N Y i=1 Mi(qi) = N Y i=1 exp(qiBi),(15) where qis the current joint configuration and Biare screw axes of the joints. The geometric Jacobian JG(q)∈B1×N⊂ G1×N 4,1is a bivector valued multivector matrix and can be found as JG=B′ 1. . . B′ N,(16) where the bivector elements can be found as B′ i= i Y j=1 Mj(qj)Bi i Y j=1 f Mj(qj).(17) Twists Vand wrenches Ware also part of the algebra and hence both can be transformed in the same manner as the geometric primitives using (14). Note that, contrary to classic matrix Lie algebra, no dual adjoint operation is needed to transform wrenches. There is, however, still a duality relationship between twists and wrenches, which can be found via multiplication with the conjugate pseudoscalar Ic=Ie0[68]. Both twists and wrenches are bivectors and the space of wrenches can be found as W ∈ span{e23,e13,e12,e01,e02,e03}.(18) The inner product of twists and wrenches V · W =−pyields a scalar, where pis the power of the motion. Similarly, the inner product of a screw axis and a wrench B·W =−τyields a torque τ, which we will use for the task-space impedance control in this article. IV. METHOD We present our closed-loop tactile ergodic coverage method in three parts: (i) surface preprocessing; (ii) tactile coverage; and (iii) robot control. The surface preprocessing computes the quantities that need to be calculated only once when the surface is captured. Tactile coverage generates the motion commands for the virtual coverage agent using the precomputed quantities from the surface preprocessing and the robot controller tracks the generated motion commands with a manipulator using impedance control. 5 A. Problem Statement We formulate a tactile ergodic controller that covers target spatial distributions on arbitrary surfaces. Similarly to HEDAC, we propagate the information encoding the coverage objective by diffusing the source term. However, we utilize the non-stationary ( ˙u= 0) diffusion equation ˙u(x, τ)=∆Mu(x, τ),(19) as it allows control over the desired smoothness [69]. Since the diffusion equation depends on time, we introduce an additional time variable τ. The diffusion time τis independent of the coverage time tused by the HEDAC algorithm and unlike HEDAC, we require an initial condition u(x,0). We set the initial condition using the source term given in Equation (1) which encodes the coverage objective at the t-th timestep of the coverage, i.e., u(x,0) = s(x, τ). Additionally, here we use ∆M, which generalizes the Laplacian for Euclidean spaces ∆ to non-Euclidean manifolds M. This operator ∆Mis also known as Laplace-Beltrami operator but for conciseness we will use the term Laplacian. Our coverage domains are curved surfaces (i.e. 2-manifolds) and we capture the underlying manifold Mas a point cloud Pcomposed of nPpoints using an RGB-D camera P:= ((xi,ci) xi∈R3,ci∈ {0,...,255}3 for i= 1, . . . , nP),(20) where xiis the position of the i-th surface point in Euclidean space and ciis the vector of RGB color intensities. We assume there is a processing pipeline (i.e., such as [20], [62], [70]) which maps the point positions and colors to the probability mass piof the spatial distribution encoding the coverage objective. Accordingly, our coverage target becomes a discrete spatial distribution p(xi) = pion the point cloud P. In order to solve (19) on irregular and discrete domains, such as point clouds, we discretize the problem in space and time. Hence, we use ui,τ to denote the value of the potential field at the i-th point at the τ-th timestep. We omit the subscript iif we refer to all points. B. Surface Preprocessing First, we compute the spatial discretization of the Laplacian ∆M. Note that there are various approaches for discretizing the Laplacian on point clouds [63]–[66]. In this work, we follow the approach presented in [66] and show a simplified version of it here, but refer the readers to the original work for more details. Using this method, the discrete Laplacian is represented by the matrix L∈RnP×nP L=M−1C,(21) where Mis the diagonal mass matrix and Cis a sparse symmetric matrix called the weak Laplacian. The entries of Mcorrespond to the Voronoi cell areas in the local tangent plane around each point of P. Similarly, the entries of C are determined by the connectivity of the points on the local tangent space and the distance between the connected points. Note that the local tangent space structure also identifies the boundary points. For a given point, the lines between the original point and its neighbors are constructed. If the angle between two consecutive lines is greater than π/2, the point is a boundary and its boundary condition is set as zero-Neumann. Next, we discretize the diffusion equation (19) in time and incorporate the discrete Laplacian L. Using the backward Euler method, we derive the implicit time-stepping equation, which remains stable for any timestep τ 1 τ(uτ−u0) = Luτ,(22) where u0and uτare column vectors containing the potential field values at the vertices of the point cloud at the initial and final times, respectively. Then, combining (21) and (22) and solving for uτwe obtain the linear system uτ= (M−τC)−1Mu0.(23) Note that solving (23) requires inverting a large sparse matrix, which might be computationally expensive depending on the size of the point cloud and requires the timestep to be set before the inversion. Alternatively, we can solve the problem in the spectral domain by projecting the original problem and reprojecting the solution back to the point cloud. This procedure generalizes using the Fourier transform for solving the diffusion equation on a rectangular domain in Rn to arbitrary manifolds. Note that the Fourier series are the eigenfunctions of the Laplacian ∆in Rn. Therefore, we can use the eigenvectors of the discrete Laplacian Lfor solving the diffusion equation on point clouds. We can write the generalized (i.e., M=I) eigenvalue problem for the Laplacian as Cϕm=λmMϕm,(24) where {λm,ϕm}are the eigenvalue/eigenvector pairs. Since Mis diagonal and Cis symmetric positive definite, by the spectral theorem, we know that the eigenvalues are real, nonnegative, and in ascending order analogous to the frequency. Therefore, we can use the first nMeigenvalue/eigenvector pairs as a low-frequency approximation of the whole spectrum. Furthermore, the eigenvectors are orthonormal with respect to the inner product defined by the mass matrix M. Accordingly, we can stack the first nMeigenvectors ϕmas column vectors to construct the matrix Φ∈RnP×nMencoding an orthonormal transformation Φ⊤MΦ=I. Then, we can transform the coordinates (shown with superscripts) from the point cloud to the spectral domain uϕ=Φ⊤Mux.(25) Note that this step is equivalent to computing the Fourier series coefficients of a target distribution in SMC. Due to the orthonormal transformation, the PDE on the point cloud becomes a system of decoupled ODEs in the spectral domain. It is well known that the solution of a first-order linear ODE ˙x(τ) = −cx(τ)is given by x(τ) = e−cτ x(0), where cis a constant and x(0) is the initial condition. Therefore, the solution of the system of ODEs in the spectral domain is given in matrix form as uϕ τ=e−λ1τ. . . e−λmτ⊤⊙uϕ 0,(26) 6 where ⊙denotes the Hadamard product. We observe from (26) that the exponential terms with larger eigenvalues (i.e., higher frequencies) will decay faster. Therefore, approximating the diffusion using the first nMcomponents introduces minimal error. Secondly, similar to the mixed norm used in SMC, the low-frequency spatial features are prioritized. Next, we transform the solution back to the point cloud to get the diffused potential field ux=Φuϕ.(27) We can combine (25), (26) and (27) into a unified spectral scheme uτ= Φ e−λ1τ. . . e−λmτ⊤⊙(Φ⊤Mu0).(28) We omit the superscripts when working on the point cloud for brevity. Note that τis the only free parameter in the diffusion computation. However, its value should be adapted according to the mean spacing between the adjacent points hon the point cloud. For that purpose, we introduce the hyperparameter α > 0and embed it into the timestep calculation τ=αh2.(29) Accordingly, we can control the diffusion behavior independently of the point cloud size. Increasing αresults in longer diffusion times and attenuates the high-frequency spatial features (see (26) for details). This corresponds to a more global coverage [69]. Conversely, decreasing αresults in shorter diffusion times, which leads to preserving the high-frequency spatial features, hence more local coverage behavior. Note that the Laplacian is determined completely by the connectivity on the local tangent space and the distance between these connected points. Therefore, it is invariant to distance preserving (i.e., isometric) transformations such as rigid body motion or deformation without stretching. Accordingly, we compute C,Mand derived quantities only once in the preprocessing step for a given surface. Recomputation is not necessary if the object stays still, moves rigidly, or the target distribution pichanges. C. Tactile Ergodic Coverage We model the actual coverage tool/sensor as a compliant virtual coverage agent shaped as a disk with radius ra. Notably, one can represent arbitrary tool/sensor footprints as a combination of disks [69]. We position our agent at the endeffector of our manipulator. Thus, for a given kinematic chain and joint configuration q, we can use the forward kinematics to compute the position of our agent as a conformal point Pa Pa=M(q)e0f M(q).(30) Since the point cloud is discrete and the agent should move continuously on the surface, we project our agent Paand its footprint to the closest local tangent space on the point cloud. 1) Local Tangent Space and Coverage Computation: Given the agent’s position Pa, we first compute the closest tangent space on the point cloud. For that, we query a K-D tree T(P) for the points xi∈ P that are within the radius raof the agent. Then, we compute the conformal embeddings Piof the neighboring Euclidean points xiusing (9). We refer to the set composed of points Pias the local neighborhood. Then, we fit a tangent space to the local neighborhood by minimizing the classical least squares objective min nN X i=1 (Pi·X∗)2,(31) where X∗is the dual representation of either a plane or a sphere and the inner product ·is a distance measure. In CGA, planes can be seen as limit cases of spheres, i.e. planes are spheres with infinite radius. This is also easy to observe by looking at Equations (12) and (13) which construct these geometric primitives. Note that fitting a local tangent sphere with the radius determined by the local curvature would always result in smaller or equal residuals than fitting a plane. It has been shown in [71] that the solution to the least squares problem given in (31) is the eigenvector corresponding to the smallest eigenvalue of the 5×5matrix bj,k = nN X i=1 wi,jwi,k,(32) where wi,k =     pi,k if k∈ {1,2,3} −1if k= 4 −1 2p2 iif k= 5. (33) Using the five components viof this eigenvector we can find the geometric primitive as X= (v0e0+v1e1+v2e2+v3e3+v4e∞)∗.(34) Note that if Xis a plane then v0= 0, otherwise Xis a sphere. Next, we want to project Pato Xby using the general subspace projection formula of CGA Ppair =(Pa∧e∞)·XX−1.(35) Here we first construct the pointpair Pa∧e∞, where e∞ corresponds to the point at infinity. Pa∧e∞is also called a flat point. Note that the projection essentially amounts to first constructing the dual line (Pa∧e∞)·Xthat passes through the point Paand is orthogonal to X, then intersecting this line with the primitive X. If Xis a sphere, then the intersection of the line and the sphere will result in two points on the sphere. If Xis a plane, it will result in another flat point, i.e. one point on the plane and one at infinity. In any case, we can retrieve the closer one to the agent position Pausing the split operation P′ a=split [Pp].(36) Here, P′ ais the projected agent position on the tangent space X. Next, we compute our agent’s footprint (i.e., instantaneous coverage) by projecting its surface to the point cloud. If the target surface was flat, all the points within the radius raof 7 our agent P′ awould be covered by the footprint. However, in the general case, both the tool and the surface can be curved and deformable. For simplicity, we assume that the surface is rigid, and it deforms the tool with a constant bending radius. We use the radius of the local tangent sphere that we computed using CGA as an approximation for the bending radius. Accordingly, we can quantify the error of the local tangent space approximation for the i-th neighbor Piby the normalized residuals eiof the least squares computation (31). We encode this approximation error into the footprint by weighting the i-th neighbor by the Gaussian kernel φ(r)using the normalized residuals ri=ei/max(e) φ(ri) = exp −ε2r2 i,(37) where the hyperparameter ε > 0controls the coverage falloff. Next, we substitute the Gaussian kernel weighted footprint into (3) to compute the coverage ct, which is then used to calculate the virtual source term stvia (1). As mentioned earlier, this virtual source term serves as the initial condition for the diffusion equation (19) at each iteration of the tactile coverage loop, i.e., u0=st. 2) Gradient of the Diffused Potential Field: We guide the coverage agent using the gradient of the diffused potential field as the acceleration command ¨ P′ a=∇uP′ a,τ ,(38) where ∇uP′ a,τ denotes the gradient of the diffused potential field at the projected agent position P′ a. However, computing the gradient on the point cloud is more involved than a regular grid or a mesh. Recall that in Section IV-C1, we already computed the projected agent position P′ a, the local neighborhood and the tangent space X∗. As the first step, we compute the tangent plane Ea,τ at P′ a, namely Ea,τ =L∗ a,⊥∧P′ a∧e∞,(39) using the line La,⊥, which is orthogonal to the surface and passes through P′ a. It is found by wedging the dual primitive Xwith P′ ato infinity with La,⊥=X∗∧P′ a∧e∞.(40) Then, we project the points Piin the local neighborhood to the tangent plane Ea,τ using (35) and (36), by setting Ea,τ as the primitive X. Next, we use the values of the potential field at the neighbor locations as the height hi=ui,τ of a second surface from the tangent plane. Then, we fit a 3-rd degree polynomial to this surface as shown by using the weighted least squares objective ˆ A= arg min Atr (Y−XA)⊤W(Y−XA),(41) with the diagonal weight matrix W W=diagφ(r1), φ(r1),...φ(rm),(42) whose entries are given by the Gaussian kernel (37). One can refer to [72] for the details. Lastly, we calculate the gradient at the projected agent’s position using the analytical gradients of the polynomial. We depict the approach visually in Figure 2. Fig. 2: Blue-red points show the value of the potential field uτ on the pointcloud Pand the yellow point is the projected agent position P′ a. We also project the agent’s neighbors Pito the tangent plane Ea,τ , shown in green. Next, we use the height function hi=ui,τ which uses the values of the potential field to lift the projected points in the normal direction of the tangent plane. We show the lifted points with large blue-red points. We fit a polynomial to this lifted surface and compute its analytical gradients at the neighbor locations ∇ui,τ , as shown with arrows in the detail view. D. Robot Control There are several aspects that the control of the physical robot needs to achieve. The first is to track the virtual coverage agent on the target surface, while keeping the end-effector normal to the surface. The second is to exert a desired force on the surface. To do so, we design a task-space impedance controller while further exploiting geometric algebra for efficiency and compactness. The control law is of the following form τ=−J⊤· W,(43) where J∈B1×N⊂G1×N 4,1is the Jacobian multivector matrix with elements corresponding to bivectors, Wis the desired task-space wrench and τare the resulting joint torques. Before composing the final control law, we will explain its components individually. 1) Surface Orientation: From Equation (40), we obtained a line La,⊥that is orthogonal to the surface that we wish to track. In [73], it was shown how the motor between conformal objects can be obtained. We use this formulation to find the motor between the target orthogonal line and the line that corresponds to the z-axis of the end-effector of the robot in its current configuration, which is found as Lee =M(q)(e0∧e3∧e∞)f M(q).(44) Then, the motor MLeeLa,⊥, which transforms Lee into La,⊥ can be found as MLeeLa,⊥=1 C(1 + La,⊥Lee),(45) where Cis a normalization constant. Note that Cdoes not simply correspond to the norm of 1 + La,⊥Lee, but requires a more involved computation. We therefore omit its exact computation here for brevity and refer readers to [73]. We can now use the motor MLeeLa,⊥in order to find a control command for the robot via the logarithmic map of motors, i.e. VLa,⊥= log MLeeLa,⊥.(46) 8 Of course, if the lines are equal, MLeeLa,⊥= 1 and consequently VLa,⊥= 0. Note that VLa,⊥is still a command in task space (we will explain how to transform it to a joint torque command once we have derived all the necessary components). Another issue is that algebraically, VLa,⊥corresponds to a twist, not a wrench. Hence, we need to transform it accordingly. From physics, we know that twists transform to wrenches via an inertial map, which we could use here as well. In the context of control, this inertia tensor is, however, a tuning parameter and does not actually correspond to a physical quantity. Thus, in order to simplify the final expression, we will use a scalar matrix valued inertia, instead of a geometric algebra inertia tensor and choose to transform the twist command to wrench command purely algebraically. As it has been shown before, this can be achieved by the conjugate pseudoscalar Ic=Ie0[68]. It follows that WLa,⊥=VLa,⊥Ic,(47) and WLa,⊥now algebraically corresponds to a wrench. 2) Target Surface Force: Since this article describes a method for tactile surface coverage, the goal of the robot control is to not simply stay in contact with the surface, but to actively exert a desired force on the surface. First of all, we denote the current measured wrench as Wm(t)and the desired wrench as Wd. Both are bivectors as defined by Equation (18). We use Wdw.r.t. end-effector in order to make it more intuitive to define. Hence, we need to transform Wm(t)to the same coordinate frame, i.e. W′ m(t) = f M(q)Wm(t)M(q).(48) In order to achieve the desired interaction force, we simply apply a standard PID controller in wrench space, i.e. WC=Kp,WWe+Ki,WZ⊤ 0 We(τ)dτ+Kd,W d dtWe(t), (49) where the wrench error is We(t) = Wd− W′ m(t),(50) where Kp,W,Ki,Wand Kd,Ware the corresponding gain matrices, and WCis the resulting control wrench. Since the desired wrench is defined in end-effector coordinates, it usually amounts to a linear force in the z-direction of the end-effector frame, i.e. Wd=fde03. Additionally, for an improved cleaning behavior one could also set a desired torque around that axis by adding τde12. The pattern of how to set this torque, however, would be subject to further investigation. 3) Task-Space Impedance Control: Recalling the control law from Equation (43), we now collect the terms from the previous subsections into a unified task-space impedance control law. We start by looking in more detail at the Jacobian J. Previously, we mentioned that we are using the current endeffector motor as the reference, hence, we require the Jacobian to be computed w.r.t. that reference. This is therefore not the geometric Jacobian that was presented in Equation (16), but a variation of it. The end-effector frame geometric Jacobian Jee Gcan be found as Jee G=Bee 1. . . Bee N,(51) where the bivector elements can be found as Bee i=f Mee i(q)BiMee i,(52) with Mee i= i Y j=N Mi(qi).(53) Hence, the relationship between JGand Jee Gcan be found as Jee G=f M(q)JGM(q).(54) The wrench in the control law is composed of the three wrenches that we defined in the previous subsections. As commonly done, we add a damping term that corresponds to the current end-effector twist and as before, we transform it to an algebraic wrench, i.e. WV=Jee G˙qe0∞.(55) With this, we now have everything in place to compose our final control law as τ=−Jee,⊤ G·KLa,⊥WLa,⊥−DVWV+WC,(56) where KLa,⊥is a stiffness and DVa damping gain. V. EXPERIMENTS Our experimental setup comprises a BotaSys SensOne 6axis force torque (F/T) sensor attached to the wrist of a 7axis Franka robot manipulator and a custom 3-D printed part attached to the F/T sensor. The custom part interfaces an Intel Realsense D415 depth camera and a sponge at its tip. We consider the sponge’s center point to be the coverage agent’s position Pa. Before the operation, we perform extrinsic calibration of the camera to combine the depth and RGB feeds from the camera and to obtain its transformation with respect to the robot joints. Additionally, we calibrate the F/T sensor to compensate for the weight of the 3-D printed part and the camera. We show the experimental setup on the left of Figure 1. A. Implementation Details The pipeline of our tactile ergodic coverage method consists of three modules: (i) surface acquisition, (ii) surface coverage and (iii) robot control. Figure 3 summarizes the information flow between the components. 1) Surface Acquisition: The surface acquisition node is responsible for collecting the point cloud and performing preprocessing operations described in Section IV-B. We use scipy1for the nearest neighbor queries and for solving the eigenproblem in (24). The matrices Cand Mcomposing the discrete Laplacian in (21) are computed with the robust_laplacian package2[66]. 2) Surface Coverage: The surface coverage node performs the computations based on the procedure given in Section IV-C. It uses the information provided by the surface acquisition node and produces the target line for the robot control node. 1https://scipy.org 2https://github.com/nmwsharp/robust-laplacians-py 9 Surface Acquisition Surface Coverage Robot Control Inner loop Outer loop Target, K-D Tree, Laplacian eigenbasis Actual coverage Line object Actual agent position Force target Fig. 3: Information flow between the three components. The pipeline is composed of an outer loop responsible for controlling the coverage progress with the feedback from the camera, whereas the inner loop compensates for the mismatch due to the robot dynamics. 3) Robot control: On a high level, the robot control can be seen as a state machine with three discrete states. The first two states are essentially two pre-recorded joint positions in which the robot is waiting for other parts of the pipeline to be completed. One of these positions corresponds to the picture-taking position, i.e., a joint position where the camera has the full object in its frame and the point cloud can be obtained. The robot is waiting in this position until the point cloud has been obtained, afterwards it changes its position to hover shortly over the object. In this second position, it is waiting for the computation of the Laplacian eigenfunctions to be completed, such that the coverage can start. The switching between those two positions is achieved using a simple joint impedance controller. The third and most important state is when robot is actually controlled to be in contact with the surface and to follow the target corresponding to the coverage agent. This behaviour is achieved using the controller that we described in Section IV-D. The relevant parameters, that were chosen empirically for the real-world experiments, are the stiffness and damping of the line tracking controller, i.e. KLa,⊥=diag(30,30,30,750,750,300) and DV= diag(10,10,10,150,150,50), as well as the gains of the wrench PID controller, i.e. Kp,W= 0.5,Ki,W= 5 and Kd,W= 0.5. The controller has been implemented using our open-source geometric algebra for robotics library gafro3that we first presented in [74]. Note that in some cases, matrixvector products of geometric algebra quantities have been used for the implementation, where the mathematical structure of the geometric product actually simplifies to this, which can be exploited for more efficient computation. B. Simulated Experiments 1) Computation Performance: In order to assess the computational performance, we investigated the two main operations of our method: (i) preprocessing by solving either the eigenproblem (24) or matrix inversion in (23) (ii) integrating the diffusion at runtime using either the spectral (28) or implicit (23) formulations. In this experiment, we used the Stanford Bunny as the reference point cloud and performed voxel filtering to set the point cloud resolution. We present the 3https://gitlab.com/gafro results for the preprocessing in Figure 4 and for the runtime in Figure 5. 3k 6k 9k 12k 15k 18k 5 0.1 2 5 1 2 5 10 2 5 25 50 100 200 implicit Fig. 4: Computational complexity of the preprocessing step for different nPand nM. Legend shows nMvalues. The time axis is logarithmic and the legend shows nMvalues. 3k 6k 9k 12k 15k 18k 0.1 2 5 1 2 5 10 2 5 100 2 25 50 100 200 implicit Fig. 5: Computational complexity of integrating the diffusion equation at runtime for different nPand nM. The time axis is logarithmic and the legend shows nMvalues. 2) Coverage Performance: We tested the coverage performance in a series of kinematic simulations. As the coverage metric, we used the normalized ergodicity over the target distribution, which compares the time-averaged statistics of agent trajectories to the target distribution εt=∥max (p−ct,0) ∥2 PnP i=1 pi .(57)