Contact Modelling and Tactile Data Processing for Robot Skins

Tactile sensing is a key enabling technology to develop complex behaviours for robots interacting with humans or the environment. This paper discusses computational aspects playing a significant role when extracting information about contact events. Considering a large-scale, capacitance-based robot skin technology we developed in the past few years, we analyse the classical Boussinesq–Cerruti’s solution and the Love’s approach for solving a distributed inverse contact problem, both from a qualitative and a computational perspective. Our contribution is the characterisation of the algorithms’ performance using a freely available dataset and data originating from surfaces provided with robot skin.


Introduction
The problem of characterising the physical interaction between robots and humans or the environment, typically using an artificial sense of touch, has received increasing attention in the literature [1][2][3]. Extensive work has been done to allow robots to obtain information about contact events [4]. Large-scale, whole-body robot skin is a key enabling technology to implement interactive robot behaviours, specifically to provide control algorithms with reliable information about contact features [5,6].
Two classes of approaches for obtaining meaningful information about contact events can be identified. The first adopts data-driven, machine learning frameworks to deal with situations where it is not trivial to model the underlying transduction principles. As pointed out in Ref. [7], data-driven approaches are appealing when modelling contact events is complex and it is difficult to model the sensor's response, specifically to take noise into account. The second focuses on appropriately modelling physical laws. Model-driven approaches are adopted whenever a (possibly simplified) model of the force distribution is available. Such models are usually based on principles of contact mechanics, and are aimed at determining closed-form solutions for contact shape reconstruction [8][9][10]. As a consequence, the two approaches originate qualitatively different results, i.e., typically a class label in the first case or a contact shape in the second.
For both data-driven and model-based approaches, three requirements must be considered: R 1 Generalisation: Run-time contact reconstruction or la belling should not assume any a priori Model of the object in contact. R 2 Efficiency: The reconstruction or classification algorithm's execution time should be predictable. Reconstructing the contact shape on the robot's surface or classifying it to inform robot behaviours require finding a trade-off between all requirements above, which may be in contrast with each other. They have been considered only to a limited extent in the literature. Table 1 reports selected data-driven and model-based approaches considering such an interplay. It indicates the worst-case computational complexity of each method using Big O notation, where c is the number of classes to discriminate from, n the size of the training set, d the data size, and f is the number of features. With the sole exception of [11], all approaches target a subset of the requirements, and do not consider them as a whole. If we restrict the analysis to the approaches aimed at reconstructing the contact shape and considering requirements R 1 , R 2 or R 3 , only the work in Ref. [7,11] appears relevant. Table 1. Classification of selected contact modelling approaches.
Kim et al. [12] c 2 Goger et al. [13] ndc Tawil et al. [14] d 2 Drimus et al. [15] ndc Decherchi et al. [16] n f Liu et al. [17] n Bhattacharjee et al. [18] ndc Ho et al. [19] d 2 + d 2 log(d) Xu et al. [20] d 2 Muscari et al. [11] d 2 Seminara et al. [7] d 2 It is necessary to better characterise the interplay between generalisation, computational efficiency and scalability when solving the problem of reconstructing the contact shape using tactile information. The contribution of this paper is a discussion about the application of foundational model-based approaches, such as the Boussinesq-Cerruti [21,22] and the Love solutions [23][24][25] to the problem of reconstructing the contact shape when large-scale, capacitance-based robot skin is used, from a generalisation, computational efficiency and scalability perspectives. A software framework with algorithms to perform specific tests is available online as open source software [26].
The paper is organised as follows. Section 2 describes the problem we want to solve using the Boussinesq-Cerruti and the Love formulations. Reconstruction algorithms are described in Section 3. A discussion about the physical plausibility of solutions is reported in Section 4. Implementation details and performance issues are described in Section 5. Conclusions follow.

Assumptions and Problem Statement
We target the ROBOSKIN technology [4,27]. A single transducer (i.e., a taxel) t i is a layered structure: the bottom layer Γ b is a positive electrode, the top layer Γ t is a ground electrode, and the mid layer is a soft elastomeric material [28], see Figure 1. In ROBOSKIN, taxels are spatially arranged in 3 cm side triangular modules, which can be connected to form large patches. Each module hosts 12 taxels. A normal force F exerted on Γ t produces variations in each taxel capacitance: where C c,i and h c,i are, respectively, the capacitance value and the elastomer thickness for t i when contact occurs, C n and h n correspond to the no contact or nominal case, 0 is the dielectric constant, r is the relative static permittivity, and A is the taxel's area. Contact induces an increase in capacitance, i.e., ∆C i > 0 corresponds to a contact event. The top layer is deformed by a pressure distribution, which is measured by taxels t i .
We are interested in solving an Inverse Elastic Problem (IEP), i.e., reconstructing the contact shape originating from the application of a force distribution [29]. This implies determining the shape of Γ t (in terms of surface tractions) using taxel readings in Equation (1). We resort to the theory of linear elasticity, and we pose the following assumptions [30]: (i) there exists a linear relation between stress and strain (deformation), and (ii) strains are infinitesimal. As discussed in Ref. [28], these assumptions are reasonable for ROBOSKIN: Stress-strain relationships are shown to be piecewise linear for different elastomeric materials. Therefore, the superposition principle can be applied, i.e., the strain resulting from a set of stresses is given by the sum of strains caused by individual stresses. Solutions to problems grounded in the theory of linear elasticity deal with elastic half-spaces, i.e., solids bounded by a plane conventionally defined by z = 0. Due to its finite thickness, ROBOSKIN cannot be treated as an open half-space. However, using the superposition principle, strains can be modelled using effective surface displacements. For a given taxel t i , if u(σ i ) is the deformation caused by the stress σ i , the corresponding effective surface displacement δ i is: where h c,i is the elastomer thickness for taxel t i . Therefore, if Q ∈ R M is a vector of surface tractions and D ∈ R N a vector of displacements, IEP consists in finding a function g such that: In our case, the function g is linear and requires the inversion of a matrix C ∈ R M×N . The use of displacements allows us to solve IEP using the Boussinesq-Cerruti or the Love formulations. Both formulations are well-known, and practical information about to solve them, including how to invert C, have been discussed elsewhere [7,29]. Here, we focus on the assumptions needed to implement algorithms considering the requirements introduced above.
Normal forces. Since ROBOSKIN adopts capacitance-based sensors, it can detect only normal forces. While this limits the type of contact events that can be detected, it allows for the simplification of IEP. Only a subset of tractions (respectively, displacements) are to be included in Q (respectively, in D), i.e., the elements of C related to tractions along the x and y axes in the elastic half-space can be explicitly set to 0. As discussed in Ref. [7], this makes C a sparse matrix, which can be efficiently inverted in O(M 2 N) [31].
Taxels layout. The solution to IEP must be independent from the underlying robot skin taxels layout, which depends on design and manufacturing considerations [32]. We discretise the IEP domain in an array of spatially distributed tractions and displacements by means of a virtual grid. Each grid cell contains a discretised value corresponding to tractions Q or displacements D in the cell's location. Although it is not explicitly required in our approach, we will consider grids in R 2 . This requires us to obtain a 2D representation of the robot's surface provided with robot skin. Methods to achieve this have been discussed in Refs. [5,33]. It is noteworthy that grid cells are the codomain of a chart whose domain is in 3D space, and in this case the manifold is the robot's surface.

On the Reconstruction of Contact Shapes from Robot Skin Measurements
In this work, we are interested in solving an IEP in which surface tractions Q represent the contact shape, and displacements D are given as a function of taxel's measurements in the form of Equation (1). First, we introduce a solution based on the Boussinesq-Cerruti model [21,22], then we point out its limitations, and finally we introduce the Love's model [23][24][25], which we adapt to our robot skin design.

Boussinesq-Cerruti's Formulation
Derivation of Influence Coefficients. The formulation of IEP according to the Boussinesq-Cerruti solution assumes the displacements in the elastomeric layer to be caused by a distribution of concentrated forces acting at discrete locations on its surface [29]. In particular, the Boussinesq-Cerruti solution for IEP is related to both concentrated normal and tangential forces F. For elastomeric materials whose Poisson's ratio is equal to ν = 0.5, which is a good approximation for the elastomer used in the employed ROBOSKIN prototypes (i.e. Ecoflex [28,34], ), the original formulas determining deformation contributions along x, y and z simplify significantly to: where: and E is the Young's modulus of the dielectric material.
Tractions vector Q consists of N = 3L components q j of L concentrated forces F l , and each concentrated force has three spatial components, respectively along x, y and z: These concentrated forces result in a vector of displacements D along x, y and z at M = 3K distinct locations: If we consider effective surface displacements Equation (2), a 3 × 3 sub-matrix of C containing influence coefficients, which relate effective displacements at location k with the forces applied at location l, is determined: wherek ranges between 3k and 3k + 2, andl ranges between 3l and 3l + 2. Influence coefficients in Equation (10) can be expressed in terms of geometric and physical properties of the elastomer, using Equations (5) and (6): where the term x kl (respectively, y kl ) is the distance between the application point of force F l and the location of the displacement δ k projected onto the x (respectively, y) axis.
Limitations of the Boussinesq-Cerruti's Solution.A brief analysis of the coefficients in Equation (10) shows that there exists a singularity in the solution for x kl = y kl = 0. Specifically, five coefficients become infinite for that choice of parameters. Geometrically, this choice corresponds to the situation in which the point of application of a force vector F l on the robot skin and the measurement's location of the displacement δ k coincide.
To overcome this singularity issue, an approximate solution has been proposed in Ref. [11] (The interested reader is referred to the referenced publication for more details. Only a conceptual overview is given here). In order to obtain it, Muscari and colleagues pretend that the employed robot skin has a thickness of h c + z 0,h c where z 0,h c is a small offset. Forces are considered to be applied to the top of the approximated robot skin, and the corresponding strain is computed using the work-force theorem and the Hooke's law. The final form for the generic displacement δ k of the approximate solution is: where Ψ(x) = x −1 (0.2431x − 0.1814) and: with A(s k ) is the area over which the concentrated force F l is exerted, i.e., in our implementation taken to be equal to the area of a grid cell s k . In the original implementation, it is taken that since z 0,h c h c , the Ψ(x) ≈ 0.25 and it is explicitly set to this value (see Figure 2 for a comparison). The solution is presented along the x axis with y = z = 0, for a force F z = 1 N, elastic modulus E = 2.1 × 10 5 Pa, robot skin thickness h n = 2 mm, and size of the grid cell dx = dy = 0.2 mm. It is noteworthy that the theoretical displacement at It is noteworthy that Muscari and colleagues consider only regular grids with equal grid cells, and therefore, assume a uniform spacing between the nodes representing points of application of force and sensing locations. Since these nodes are placed perfectly one over another, they apply the approximate solution when computing influence coefficients for nodes that correspond to the same point in space. In other cases, it is obvious that either one of x kl or y kl is non-zero, and therefore, the solutions are not singular. It logically follows that the criterion the authors apply for choosing between the direct Boussinesq-Cerruti solutions and the approximated solution is discrete in nature. In order to apply the Boussinesq-Cerruti solutions to any combination of (not necessarily regular) grids, with arbitrary cell sizes, a continuous criterion must be employed instead. To this end, one must consider the nature of the functions constituting the Boussinesq-Cerruti solutions. The singularity occurs with z = 0, and therefore, the problem is planar in nature, i.e., solutions are functions of (x kl , y kl ). If these functions could be shown to be concave (i.e., their Hessian matrices were positive semi-definite), then it would be trivial to choose a solution ( Figure 3). Unfortunately, this cannot be shown trivially. If we consider the expression of the influence coefficient c 3k+2,3l+2 in Equation (10) for a generic displacement δ z j and force F z k , we express it adopting a polar coordinate system r 2 , φ , where r 2 = x 2 kl + y 2 kl and φ = atan2(y kl , x kl ), and we try to find critical points by calculating: it can be shown by applying the Sturm's theorem to the resulting polynomial that there exists a critical point in the interval r ∈ (0, ∞). However, determining the nature of the critical point would require us to solve a quartic polynomial and possibly further steps to examine growing orders of derivatives. This approach was deemed unsatisfactory, but no better continuous criterion can be proposed at this time. For this reason, the aforementioned criterion is what the current implementation employs, i.e., both the original Boussinesq-Cerruti and the approximated solutions are computed, and then the one with the lower absolute value is used for the influence coefficient.

Love's Formulation
Derivation of Influence Coefficients. The formulation of the elastic problem for the Love's solution considers displacements δ k in the elastomeric material to be determined by a distribution of normal pressures acting over rectangular areas on the robot's surface [24]. Such areas, or grid cells, have sizes 2a and 2b ( Figure 4). For the sake of notation, let us assume that a point on the boundary of the half-space (i.e., the robot skin's surface Γ t ) has Cartesian coordinates x , y , 0 and a point within the elastomer has Cartesian coordinates x, y, z with z > 0. Furthermore, let r be an auxiliary variable such that r 2 = ∆x 2 + ∆y 2 + z 2 , where ∆x = x − x and ∆y = y − y. Therefore, solutions follow the notion of spatial derivatives of elastic potential functions introduced by Boussinesq, whereas deformations are computed as: where χ and V are two functions describing the elastic potentials, defined respectively as: and: where p is the pressure acting on the grid cell. The term χ in Equation (18) is called Boussinesq's 3D elastic logarithmic potential, while the term V in Equation (19) is the Newtonian potential of the surface distribution. As described in Ref. [25], integrals in Equations (18) and (19) can be evaluated to obtain closed-form solutions for the deformations in an elastic half-space excited by a uniform normal pressure applied over a rectangular grid cell, as follows: where: In Equations (23)-(30), the parameter j determines which sign is used, i.e., the upper (respectively, lower) sign corresponds to j = 1 (respectively, j = 2), and p is the value of the (uniform) pressure acting on the cell.
It is noteworthy that deformations are computed as linear functions of the pressures exerted on the robot skin and therefore, the model is still linear, which allows us to take advantage of the superposition principle. In the Love's formulation, the tractions vector Q consists of N normal components of a pressure distribution: These pressures result in a distribution of displacements D along the x, y and z directions at K locations in the same form of Equation (9). Therefore, the 3 × 1 sub-matrix of C contains influence coefficients which relate displacements at location k with pressures applied at location n. According to Equation (2), these can be expressed as: where: ln ∆y + r 20 ∆y + r 10 z=h c y =b and terms J 1 , J 2 , K 1 , K 2 , L 1 , L 2 are characterized by Equations (23)- (25).  Limitations of the Love's Solution. As it can be noticed in Figure 6, if the size of the cell is decreased, the Love's solution smooths out as well and becomes similar to the Boussinesq-Cerruti's solution, without the singularity. This suggests that Love's formulation assures a certain degree of complexity required to model highly complex contact situations, specifically when the size of the grid cells is reduced. It should also be noted that equations in the Love's formulation do not contain real singularities, which is a fundamental feature to adopt the model with real-world robot skins. As a matter of fact, there are expressions becoming singular in specific cases, but it can be shown that in all of them, they are multiplied by a coefficient which tends to 0 faster than the expression becomes singular. However, special care must be taken when implementing those equations, especially when dealing with finite precision arithmetic.

Comparison between Boussinesq's and Love's Solutions
A qualitative comparison between the two solutions presented in Section 3.1 and Section 3.2 can be done by comparing deformations resulting from similar load conditions. In order to generate similar load conditions, we shall only consider loads normal to the surface. It is noteworthy that this is precisely the case with the robot skin technology we use, and due also to the fact that the Love's problem deals with normal pressures only. Since the loads considered by the Boussinesq-Cerruti's and Love's formulations are different in nature (i.e., forces and pressures, respectively), we shall approximate the force acting in the Boussinesq-Cerruti model as: where p is the normal pressure acting in the Love's formulation, and 2a and 2b are the sizes of the grid cell (where the load is applied). Furthermore, for the approximate solution to the Boussinesq-Cerruti's model, we assume that the discretization of the forces space is done with elements such that the area over which the force is exerted is equal to that of the pressure's grid cell in the Love's model: A sample result of the comparison can be seen in Figure 6. The two models produce results within the same magnitude order, i.e., less than 2 × 10 −4 , but nevertheless there is a significant difference in the amplitude of the deflection, the largest difference occurring at x = y = 0. We can try to qualitatively compare this difference by computing the approximate Boussinesq-Cerruti's solution for a fixed force value F z and varying the area of the grid cell A(s), and plotting it with a Love's solution for grid of size 2a = 2b = A(s), and pressure value p = F z /(2a × 2b), evaluated at x = y = 0. The discrepancy in the two solutions is contrary to what one would expect. Intuitively, the Boussinesq-Cerruti's solution (i.e., concentrated force) should result in an indentation with narrower, sharper slopes than the solution corresponding to the Love's problem. On the other hand, the indentation in the point of application of the force would be expected to be deeper in the Boussinesq-Cerruti's solution, since it originates from a concentrated force, but this is not the case.

Plausibility of Solutions
By analysing the contact situation, one can derive a set of constraints on the model involved in the inverse elastic problem, which define the set of feasible solutions for tractions Q. Usually, these constraints originate from the likelihood of observing some measurements due to the laws of physics governing the contact event. It may occur that the feasible solution set does not overlap with the theoretical solutions obtained through the inversion of C. This means that there might not exist a physically feasible solution (in terms of model parameters) which would produce exactly the same data if a solution to the forward elastic problem was computed with the model's parameters taken to be such exact solutions. This phenomenon is readily explained by taking into consideration the presence of noise in the data, which is the input to the inverse elastic problem. Input data (e.g., taxel measurements) are affected by various sources of noise, most of which are unrelated to the nature of the elastomer employed in the robot skin, e.g., electric charge fluctuations on the capacitance-based taxels, limited processing precision of CDC units, or even the limited precision of computer arithmetic, just to name a few. An exact solution to the inverse problem (if it exists) corresponds to modeling all these imperfections, which is usually not desirable. One can therefore, sacrifice the least-squares guarantee of a pseudo-inverse based solution for meeting physicality constraints on the tractions set.
Following the ideas presented in Ref. [7], one can assume that under common contact conditions, normal tractions (forces or pressures) can act only in a compressive way, i.e., q z j > 0 for every j. Moreover, we assume that tangential tractions are generated exclusively by friction [35]. If we further constrain ourselves to single-contact cases only, then tangential tractions can be shown to be proportional to normal tractions.

Compressive Normal Tractions with the Fourier-Motzkin Elimination
There are many ways to constrain normal tractions to be compressive. Let us assume that there are only normal tractions acting on the robot's body. As a consequence, the relationship between displacements and tractions becomes: with a general solution in the form: where z ∈ R N such that z = [z 1 , . . . , z N ] T is an arbitrary vector. As mentioned above, constraining tractions to be non-negative can be expressed as the following system of inequalities: If we restrict the domain in which to search for the non-negative solutions to the set of exact solutions given by Equation (37), the problem of obtaining the non-negative constraints on the elements of Q is equivalent to that of finding the components of the z vector satisfying a system of linear inequalities: It can be proved that a solution in closed form to this set of linear inequalities exists, and it can be obtained using a mathematical tool introduced by Fourier in Ref. [36], later rediscovered by Lloyd L. Dines and Theodore S. Motzkin, and came to be known as the Fourier-Motzkin elimination.
The process to obtain a solution is similar to that of Gauss elimination, and consists in transforming a set of linear inequalities to an equivalent set, which is expressed without a subset of the original variables. Without loss of generality and for the sake of conciseness, let us consider an arbitrary system of linear inequalities: where the values of real coefficients a ji are not constrained, i.e., they can be positive, negative or equal to zero. Then, we can eliminate a variable x e by transforming the original inequalities so that on the left-hand side they contain only x e , and grouping them into three classes, depending on their direction (which corresponds to the sign of coefficients a je ): 1.
where b j = b j /a je and a ji = a ji /a je . The original system of inequalities is then equivalent to: In the worst case (with an equal number of inequalities in the first and second group), the number of inequalities in a system created by eliminating one variable is n 2 /4. Running p successive elimination steps will result in at most 4(n/4) 2p inequalities. Many of these inequalities are redundant and the required number of inequalities can be shown to grow as a single exponential [37].
Considering again the inverse elastic problem, without enough insight, the Fourier-Motzkin elimination might seem promising, since it could be performed as an offline step, with back-substitution being the only online part.

Remark 1.
The complexity of the Fourier-Motzkin algorithm must be analyzed. Since the displacement vector D is only known at run-time (as a consequence of a contact event detected by the robot skin), no detection of superfluous inequalities (as mentioned above) can be done offline. Therefore, the worst-case size of the final system of inequalities obtained after eliminating all variables is: which, for a moderately sized problem of reconstructing a distribution of 100 tractions (e.g., a 10 × 10 grid) corresponds to approximately 1.6 × 10 280 inequalities. Without elaborating the details, encoding this in an actual robot software architecture would require approximately 5.8 × 10 273 gigabytes of memory under double precision.
From this comparison, it should be clear that the Fourier-Motzkin elimination cannot be adopted in practice.

Compressive Normal Tractions with Non-negative Least Squares
The problem of finding non-negative solutions to systems of linear equations is a well studied problem in mathematical optimization, driven mostly by engineering applications [38]. As it turns out, quite often physically feasible solutions are constrained to the non-negative only domain.
The first algorithm to find a least-squares solution in the non-negative domain is due to Lawson and Hanson [39]. This algorithm belongs to the family of active set algorithms. It tries to iteratively optimise a subset of system variables whose values can change (i.e., the free set) so that the final solution satisfies the non-negativity constraint with a minimal least square-error. However, the performance of the algorithm proposed by Lawson and Hanson is hindered by the expensive computation of a matrix inversion related to the free subproblem at each refinement step of the solution. Moreover, this method belongs to the single pivoting group, meaning that at each iteration, only one variable is chosen to be moved between the free and the active sets, which increases the number of iterations required to reach a good solution. Since then, many other algorithms to compute the non-negative least-squares solution to a system of linear equations have been proposed, usually offering significant performance improvement over the Lawson-Hanson algorithm. Most of the algorithms have been generalised to realise box-type constraints, such that each variable x i is subject to l i ≤ x i ≤ u i constraints. A partial survey of these algorithms is given in Ref. [38]. The majority of the modern solvers belong to the group of block-pivoting, active set algorithms. The main difference with respect to the original algorithm by Lawson and Hanson is that, in each iteration, the algorithms try to move more than one variable between the free and the active sets. In general, the inverse elastic problem is dense in its nature, meaning that (almost) all elements of the influence coefficients matrix are non-zero. Nevertheless, Mikael Adlers describes an efficient and robust algorithm for solving box-constrained linear problems which remain highly efficient even for dense problems [40]. This algorithm, known as BLOCK3 block pivoting algorithm, has been chosen here to solve the inverse elastic problem with non-negativity constraints on the tractions, specifically the implementation described in Ref. [41].

Remark 2.
It is noteworthy that any algorithms implementing non-negative least squares solutions must be performed on-line, i.e., when the displacement vector D in the inverse elastic problem is known.
Although the BLOCK3 algorithm is proved to converge in a finite number of iterations, the complexity of the involved operations (and therefore, the problem of putting real-time constraints on them) may prove a limitation in the case where hard real-time constraints are required.

Implementation Details
The approach presented in this paper has been developed as an open source software. A video showing online contact shape reconstruction is available as supplementary material to this paper and it is available at this link http://www.mdpi.com/1424-8220/19/4/814/s1. The current implementation is based on C++ and is built on top of the Skinware framework [6,42,43]. In particular, we highlight the following features: • both Boussinesq-Cerruti's and Love's formulations have been implemented; • the algorithms for tractions reconstruction and surface deflections have been designed considering a clear division of the offline and online parts; • the implemented algorithms allow for the reconstruction of forces using displacements, pressures using displacements, displacements using forces, and displacements using pressures; • for the inverse problem, variants of the algorithms exist, which enforce the non-negativity constraints on tractions.
As it can be observed from the derivation of the influence coefficients (both for Boussinesq-Cerruti's and Love's approaches) above, the elements of C depend only on the geometry of the problem. This fact can be exploited to divide the solution of the reconstruction problem into an offline and online part, which is fundamental for practical reasons from a computational standpoint. When the straight least-squares solution via the pseudo-inverse is used, both the original C matrix and its pseudo-inverse C + can be computed before tactile data processing starts. Otherwise, only the C matrix can be computed in advance.
This observation can be used to greatly improve the overall performance, since the only step to perform online would become a matrix-vector multiplication, which is a relatively lightweight operation. In order to compare the time complexity of both the Boussinesq-Cerruti's and Love's solutions and to assess the gain obtained by offline computation, the C matrix has been computed for different grid sizes. Figure 7 shows computing times for an increasing number of grid cells. Influence coefficients have been computed for two identical instances of regular square grids. Only normal displacements and tractions have been considered. Therefore, the size of the C matrix is n × n, where n is the number of grid cells. As expected, Figure 7 shows a quadratic time complexity for increasing grid sizes. However, there is a large discrepancy in the growth rate for the Boussinesq-Cerruti's and Love's models. Qualitatively, the time required to compute influence coefficients by the Love's model is tenfold the time for the Boussinesq-Cerruti's model in an equivalent grid. This is possibly due to the high complexity of the formulas for computing Love's influence coefficients (Equation (25)) compared to the Boussinesq-Cerruti's solution (Equation (6)). With such an increase in the computation time, the capability of performing the majority of calculations offline becomes crucial.
Two remarks should be made. On the one hand, it is noteworthy that in case non-negativeness constraints for the reconstruction of tractions are considered, only the C matrix can be pre-computed, whereas its pseudo-inverse cannot. On the other hand, this division allows a number of computations to be performed off-line, which is a great benefit if the tractions reconstruction were to be performed in real-time. The software architecture could then enable solving the elastic problem with an adapting resampling resolution, querying external, off-site computation nodes for the required matrices. This means that resource intensive computations could be moved to a dedicated unit, possibly endowed with dedicated hardware for maximum efficiency, e.g., a GPU, whereas the use of such hardware would be otherwise prohibitive on board the robot, e.g., for reasons of power efficiency or space occupation.

Experimental Setup
Experimental data used to validate the presented models are taken from the dataset used in Ref. [11]. Data have been recorded using six triangular ROBOSKIN modules combined together to form a hexagonal patch, shown in Figure 8 on the right hand side. A three-way Cartesian robot positioner from Thorlabs Inc. allows the patch to horizontally translate along the x and y axes, and to rotate it along the vertical axis. A force of up to 3 N can be applied using a linear actuator mounted along the vertical axis. The tip of the linear actuator, which is provided with a load cell to measure the peak force on the area below, can be mechanically coupled with indenters of various shapes. The experimental setup is shown in Figure 9.  Four parameters have been varied during the experimental campaign, namely the location of the pressure exerted by the rigid indenter (i.e., over a given taxel, Figure 10a, between two taxels, Figure 10b, and between three taxels, Figure 10c), the indenter's shape (i.e., a 12 mm diameter half-sphere, a 6 mm diameter half-sphere, a 12 mm diameter cylinder, a 3 mm diameter cylinder), the value of the exerted pressure (i.e., 0.2 N, 0. 5 N, 1 N, 1.8 N, 2.5 N, and 3 N), and the duration of the contact phase (from 3 sec to 7 sec). Each trial involving any combinations of these parameters has been repeated 50 times for statistical significance.

Performance of Love's Solution Compared to Boussinesq-Cerruti's Solution
In all the experiments, Boussinesq-Cerruti's models have been computed using an approximation for the function Ψ by a constant number, i.e., Ψ = 0.25, as suggested in Ref. [11]. As a reference for further considerations, let us consider the experiment in which the robot skin patch is indented using a large sphere with a force of 1.8 N, which corresponds to Figure 11a. The resolution of the grid over which the deflections are measured is about 2 mm, see Figure 11c. A qualitative comparison of results obtained for the inverse elastic problem using the Boussinesq-Cerruti's and Love's models is shown in Figure 12. The reconstruction grid is squared, i.e., the number of tractions to reconstruct is equal to the number of available measurement data points. The two results are surprisingly similar, if we ignore the scale, and only minor differences can be noticed. This is most probably due to the fine discretisation of the grid. Each pressure acts on a 2 × 2 mm 2 size cell, which is an oversampling with respect to the actual ROBOSKIN spatial density (i.e., 2 taxels per cm 2 ).  To compare the reconstruction quality for the two solutions, surface deflections have been reconstructed: 1. determining the tractions by solving a square inverse elastic problem in which the grids have the exact same geometry and are placed one over another; 2.
reconstructing surface deflections for varying resolutions of the reconstructed displacements grid, namely (i) a grid with cells shaped as squares with a side length equal to 0.5 mm, (ii) a grid being an exact copy of the input sensor data grid, and (iii) a grid with square shaped cells with side length equal to 3 mm.
Results are shown in Figure 11 for the Boussinesq-Cerruti's solution and in Figure 13 for the Love's solution.
Several comments can be made. The first comment is that both in the Boussinesq-Cerruti's and Love's cases, reconstructing surface deflections with respect to a grid that is an exact copy of the sensory input grid, results in surface deflections which are exactly equal to the sensor input. Obviously enough, this is to be expected. In both cases, the influence coefficient matrix for the inverse problem C I is square. Moreover, the coefficient matrix describing the forward problem C F is square as well, and most importantly those two matrices are equal. Therefore, in this particular case, reconstructing effective surface displacements D corresponds to: This effect is evident when one compares Figure 11a,c, and Figure 13a,c, for Boussinesq-Cerruti's and Love's solutions, respectively.
Then, it must be remarked that Boussinesq-Cerruti's solutions perform poorly in reconstructing surface displacements for grids characterised by a resolution different from the input's, as shown in Figure 11b,d. On the one hand, it is apparent that Boussinesq-Cerruti's solutions perform well for problems in which the influence coefficient matrix is square, but once this property is lost, undesired artefacts appear. Therefore, we conclude that Boussinesq-Cerruti's solutions are not suitable for a stable resampling in the displacements field. On the other hand, if we consider the performance and stability of Love's solutions in reconstructing surface displacements for grids with a different resolution than the input, i.e., resampling the displacements field, shown in Figure 13b,d, one immediately notes the consistency of the reconstruction. For an increased resolution of the reconstruction, as shown in Figure 13b, Love's solutions yield a smooth interpolation of sensor data, which is close to what one would expect to happen in actual robot skin. If the Love's solution is used for downsampling the surface deflection field, as shown in Figure 13d, the result is a deflection field which is still representative of the contact event. It is noteworthy that the difference between the height of the main peak in the input and those of down-sampled cases is of the order of one twentieth of a millimetre. The difference between the two models is even more apparent if a downsampling of the tractions field is considered, as shown in Figure 14. In the Boussinesq-Cerruti's case, if the tractions are reconstructed on a regular grid with 3 mm grid size square cells, as shown in Figure 14a, the result is hard to interpret at best. In addition, the reconstructed surface deflections, as shown in Figure 14c, do not resemble in any way the raw input sensor data shown in Figure 11a. On the contrary, Love's solution performs better in both tractions reconstruction and surface deflections, as exemplified in Figures 14b,d. In particular, the similarity between Figures 14d and 11d must be noted. The latter Figure corresponds to reconstructing the displacements field at 3 mm, without downsampling. It shows that about the same amount of information or precision is lost regardless whether the tractions grid or the reconstructed surface deflections grid is down-sampled. This argument is backed by the dissimilarity of Figures 11c and 14d, which carry a different amount of information regardless of the fact that the reconstructed displacement grids are (approximately) of the same resolution. On the one hand, the above observations can be exploited when using the algorithms for an actual reconstruction to cut down the overall computational cost. If the aim is to down-sample the reconstructed displacements field, instead of solving an exact inverse elastic problem and then a smaller forward elastic problem, one can solve both problems at a smaller scale, since about the same amount of information is preserved, and therefore, saving precious computational time. On the other hand, increasing the actual resolution of the reconstructed tactile image by resampling the tractions field has not been explored in this work. It is noteworthy that there exists an infinity of solutions and not well-informed principles have been identified to single one out.
As shown in Section 4.2, while solving an inverse elastic problem, one may restrict the set of traction components to be non-negative, by using one of the non-negative least squares algorithms available in the literature. Figure 15 shows examples of such results. Specifically, it is noteworthy that the solutions obtained by imposing the non-negativity constraint are not very dissimilar to free solutions, for example those shown in Figure 12a,b. Furthermore, the reconstructed surface deflections resemble sensor measurements from Figure 11a with a high degree of fidelity. The already high fidelity degree of free Love's solutions, as exemplified in Figure 13, arises a question of usefulness of non-negativity constraints in the first place. The constraint does improve the solution by eliminating pathological components. However, it makes meeting hard real-time requirements harder, due to the complexity of the involved algorithms, which can be only run online. Furthermore, the online phase, which solves a system of linear equations with non-negativity constraints, is expected to be slower than in the case of a free solution, which involves merely a single matrix-vector multiplication operation. As shown in Section 3.1, the Ψ function in the approximate Boussinesq-Cerruti's solution can be approximated by a fixed number, i.e., Ψ(x) ≈ 0.25. The authors of Ref. [11] have explicitly set Ψ to this fixed value. To explore the impact of this choice, traction reconstruction with the exact calculation of the Ψ function has been performed, using both the free and the non-negativity-constrained solutions. The result is shown in Figure 16. The free solution, shown in Figure 16a, still shows a peak in the contact location, but nevertheless contains a very large amount of noise, the interpretation of which is difficult. Obviously enough, for the reasons backed by Equation (40), the reconstruction of surface displacements corresponds exactly to input data. The non-negativity-constrained solution presents an interesting filtering behaviour. Forces are purported to be exhibited only in the location of highest surface deflections. The surface deflection reconstruction is not characterised by the fidelity with respect to input data as the Love's solution. However, it is possible that such reconstruction (Ψ computed using the exact formula along with non-negativity constraints) effectively filters noise from the solution with respect to input data. Possibly, the areas with the largest discrepancies between said reconstruction and input data are those where the strongest noise is prevalent. Nevertheless, without a detailed analysis of noise sources and their nature, this statement cannot be defended trivially.

Discussions
In this section limitations and generalization aspects of the proposed software architecture are discussed.

Generalization with Respect to Other Tactile Sensors
The presented software architecture follows the model-driven approach and it receives as input a deformation (in mm) computed using the transducer model. In this work we used, as transducer, the ROBOSKIN sensor that is based on capacitive transduction as highlighted in section 2. The capacitive sensor provide a variation in capacitance which is proportional to the strain of the dielectric layer and therefore, we used the constitutive equation of a capacitor to get strain information from the output of the sensor. In this respect the proposed software architecture can be used with any transducer that could provide deformation information for example [44] and [45].

3D Surfaces
The proposed architecture considers the sensor integrated on a flat surface, this is of course a limitation since in large-scale robotic skin all the taxels can be located on the robot body so the shape of the actual sensing surface can be in 3D space. This work has been limited to 2D surface since a planar force distribution can be easily generated as a ground truth in order to benchmark the proposed computational models; furthermore, with respect to the assessment real-time reconstruction capabilities of the proposed software architecture, the performance could be considered a benchmark also for non/planar geometries.

Shear Forces
In all the presented experiments the forces have been applied normal to the sensor surface and no other force directions have been taken in account. In this respect it is important to point out that the Bousinnesq and Love's approaches take in account the tangential components of the forces, but ROBOSKIN technology, based on capacitive transduction provides information related only to strain along the normal direction therefore, the model doesn't use that information for the reconstruction of the forces. If other transducers, able to provide this missing information, are used the software architecture can provide reconstruction of shear stress.

Conclusions
This paper introduces an approach to extract information about contact events based on large-scale tactile sensors. The approach is based on a closed-form algorithm for the reconstruction of the contact shape, which is backed by a physical model of a large-scale, capacitance-based robot skin technology we have developed in the past few years, namely ROBOSKIN. The classical solution that can be obtained by adopting the Boussinesq-Cerruti's model is compared (both from a qualitative and a computational perspective) to the Love's approach for solving a distributed inverse contact problem, which has been adapted to our case. The paper elaborates on two aspects: the first is a proposal for a general-purpose algorithm for the reconstruction of deformation and force distributions in case of capacitance-based robot skins; the second is related to the characterisation of its real-time performance, which can be tuned according to available computational resources. Experiments on robot skin patches have been performed to provide a quantitative analysis of results.