{'fraction_non_alphanumeric': 0.03202880518795008, 'fraction_numerical': 0.011310937676733402, 'mean_word_length': 4.671656153105955, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theory and Design of Super-resolution Haptic Skins

Haptic feedback is important to make robots more dexterous and effective in unstructured environments. High-resolution haptic sensors are still not widely available, and their application is often bound by the resolution-robustness dilemma. A route towards high-resolution and robust skin embeds a few sensor units (taxels) into a flexible surface material and uses signal processing to achieve sensing with super-resolution accuracy. We propose a theory for geometric super-resolution to guide the development of haptic sensors of this kind and link it to machine learning techniques for signal processing. This theory is based on sensor isolines and allows us to predict force sensitivity and accuracy in contact position and force magnitude as a spatial quantity. We evaluate the influence of different factors, such as elastic properties of the material, structure design, and transduction methods, using finite element simulations and by implementing real sensors. We empirically determine sensor isolines and validate the theory in two custom-built sensors with barometric units for 1D and 2D measurement surfaces. Using machine learning methods for the inference of contact information, our sensors obtain an unparalleled average super-resolution factor of over 100 and 1200, respectively. Our theory can guide future haptic sensor designs and inform various design choices.

Haptic feedback is important to make robots more dexterous and effective in unstructured environments. High-resolution haptic sensors are still not widely available, and their application is often bound by the resolution-robustness dilemma. A route towards high-resolution and robust skin embeds a few sensor units (taxels) into a flexible surface material and uses signal processing to achieve sensing with super-resolution accuracy. We propose a theory for geometric super-resolution to guide the development of haptic sensors of this kind and link it to machine learning techniques for signal processing. This theory is based on sensor isolines and allows us to predict force sensitivity and accuracy in contact position and force magnitude as a spatial quantity. We evaluate the influence of different factors, such as elastic properties of the material, structure design, and transduction methods, using finite element simulations and by implementing real sensors. We empirically determine sensor isolines and validate the theory in two custom-built sensors with barometric units for 1D and 2D measurement surfaces. Using machine learning methods for the inference of contact information, our sensors obtain an unparalleled average super-resolution factor of over 100 and 1200, respectively. Our theory can guide future haptic sensor designs and inform various design choices.

Introduction
Haptic sensors are indispensable in robotic applications to enable robots to perceive when, where, and how their bodies are interactively contacting other things (1)(2)(3)(4)(5)(6). A common theme for surface haptic sensors is to integrate many small sensing units forming a type of grid along a flat or curved surface. Each sensing unit, named taxel, is responsible for sensing interactions near its location. For typical applications, a resolution is desirable that would imply numerous taxels (4,5,(7)(8)(9)(10)(11). This is true both for small surface sensing, e.g., at the fingertips, and for large surfaces, e.g., around limbs. For fingertips, very focused areas need a high density of taxels with a fine size to perceive high-resolution haptic information, which is similar to touch screens. For limbs, even though the haptic information is coarsely needed referring to the density distribution of human mechanoreceptors (12), but the sensing areas are comparably large which also need many taxels. Technical challenges arise concerning the physical size of the taxels as well as growing manufacturing and wiring costs (13).
Naturally, efforts have been made to reduce the number of taxels (14)(15)(16)(17)(18)(19)(20). The reduction is generally enabled by the fact that a taxel can monitor not only a tiny area but an extended patch on the surface and that we are interested in a subset of all possible stimulation patterns, for instance, a few touch-points. The particular material properties lead to a characteristic spread of contact information to the sensing taxels. Different physical effects can be used, such as geometric and mechanical properties (15)(16)(17), electrical resistance (18,19), magnetic flux (20), thermo-and fluid dynamics (21), and so forth. The central idea is to solve the inverse problem of inferring haptic information from a few sensors, effectively creating high-resolution virtual taxels. This is also referred to as super-resolution sensing (22).
Super-resolution sensing was early explored in different areas, such as geostatistics (temperature and precipitation prediction) (21,23), imaging systems (super-resolution microscopy) (24,25), and material releases detection (aerosol/chemical plume release) (26,27). They typically rely on the continuity and neighboring effects in the transmission medium and on solving inverse problems (signal processing) to reconstruct spatial information from a few sensors in a sparse configuration. There are also several haptic sensor designs that achieve super-resolution sensing.
HapDef (17) makes use of the spreading behavior of mechanical deformation by attaching a few strain gauges on a large robotic limb shell and developing machine learning (ML) methods to achieve a 78-fold super-resolution (SR). Similarly, Piacenza et al. (28) mold a few barometers inside silicone and employs ML to achieve a 57-fold SR. Lepora et al. (15) customizes a sensor with several capacitive taxels and applies Bayesian optimization to achieve a 35-fold SR. Yan et al. (20) uses Hall sensors in sparse configuration to detect the deformation of a magnetized flexible film and develops simple ML models to achieve a 30-fold SR. These sensors rely on the empirical knowledge that overlapping multiple taxels' perception fields triggers SR behavior. However, a sounded quantitative theory is still missing in solving the inverse problem and providing a basis for a fast back-projection-like approach towards system-level sensor designs. The theory would imply a baseline for evaluating signal-processing methods from a model-centric perspective (optimizing ML models based on whatever data collected), and offer guidelines of  Fig. 1. Introduction of the model. A-I depicts a 1D sensor layout with sparsely placed real taxels and highly resolved virtual taxels. A-II shows the orange taxel value isolines (TVIs) that are the force F needed to elicit the same taxel-value S depending on the distance d from the sensor unit. The intersection of TVIs of two sensors indicates a hypothetical contact point P i with a force strength F i . The measurement noise σ S in shaded blue leads to uncertainty of the position σ P i and force σ F i . B shows a validation of our theory on a 2D sensor that embeds 25 barometers in an elastomer (EcoFlex 00-30) with a grid layout (5 × 5) to realize the super-resolution functionality. As a demonstration, a 4 mm red spherical indenter contacts the square sensor surface (34 mm × 34 mm) following the "Minerva" pattern (diameter 26 mm) with different depths (1.2 mm, 2.4 mm, 3.6 mm). The sensor can resolve the contact pattern in different resolution depending on indentation depths. data collection for data-driven sensors from a data-centric perspective (choosing which data to collect given ML models).
How much "super-resolution" can we obtain? What is the expected sensitivity of the sensing system? Which surface material and which transduction method is suitable? Why and how can machine learning help the inference? This paper aims to provide answers to these questions from a theoretical perspective as well as from a practical point of view by applying the method to three common sensor types and evaluating two custom-built sensor designs. We introduce a theory that allows us to infer spatially resolved properties, such as accuracy and sensitivity for single and double contact, based on a single characteristic of the material-taxel interplay: taxel value isolines (TVIs).
The model we use is a sensor comprising discrete taxels sparsely distributed in a continuous transmission medium underneath a sensing surface as shown in Fig. 1. For sensing devices that can be described by this model, we derive the single-contact accuracy for position and force magnitude inference based on the taxel value isolines; determine the minimal force profile for which localization is possible; provide conditions for simultaneous contacts to be distinguishable; evaluate the influence of different factors, such as elastic properties of the material, structure design, and transduction methods as well as contacting object shapes; and design a one-dimensional (1D) and a two-dimensional (2D) barometer-based sensors to validate our proposed super-resolution theory.

The Model
We consider a class of haptic sensing devices intended for measuring force interactions on an extended surface that have an elastic transmission medium covering or embedding physical sensor units (taxels). A one-dimensional model is shown in Fig. 1 A-I. A single taxel value s is a function of the applied contact force strength F , and the displacement d between the contact center and the taxel center: where S is the taxel measurement noise with a constant standard deviation σ S : To relate the sensor measurement noise S to force values F , we assume a linear relationship. Formally, the sensor responses to an applied force at distance 0 should be f (F, 0) = c · F and to simplify the notation we set c = 1. However, our analysis can be easily adapted for different values of c and for nonlinear monotonous relationships. A detailed analysis of c for a sensor type used in this paper is provided in Sec. D, which supports our linearity assumption. One of our main contributions is to introduce the taxel value isolines as an important characteristic function of the system.
Definition 1 Taxel value isolines (TVIs) are a family of curves where the mean taxel output (Eq. 1) has a constant value S.
The TVI for the model system and the effect of measurement noise are shown in Fig. 1 A-II. Intuitively, the isolines quantify how much force is needed along the surface to yield the same sensor value. To activate a taxel with a particular value, the required force strength is smaller when the contact location is closer to the taxel. Based on these isolines, we can derive the accuracy distributions of contact localization, force quantification, and the distribution of sensitivity over the sensing surface. Our method is based on the conditions for unique nonlinear triangulation and the geometry of intersecting isolines with their uncertainty bands.

Super-resolution in 1D
When can a single contact point be localized at super-resolution, meaning much more accurate than the distance between taxels? Intuitively, this is possible when two or more taxels measure nonzero responses to the contact force and the activation pattern of the taxels is unique for this location. More formally, this condition can be analyzed with the TVIs. The minimal setting of two taxels located at a distance of D is shown in Fig. 2.
A particular contact event causes a sensor reading in both sensors S-1 and S-2. The TVI corresponding to a sensor reading of S-1 relates the position to the force of the potential contact point. Only if the TVIs from both sensors intersect, the contact position can be localized, up to some uncertainty introduced by the measurement noise σ S (Eq. 1). With this insight, we can derive the minimal force sensitivity F S that allows for super-resolution localization, i.e. where at least two isolines intersect. We can find this minimal force, where the isolines corresponding to the smallest taxel sensitivity S min intersects with TVIs from neighboring taxels. The area where F ≥ F S is shaded in green in Fig. 2 A-I.
As we assume an additive measurement noise, the uncertainty is also added to the isolines. What does this mean for the ability to localize the contact position and to infer force magnitude? For the case of parabolic TVIs, we have to consider different possible scenarios, as shown in Fig. 2 A-II. For contact locations between the two taxels, marked with p 1 , p 2 , the standard deviation for the force estimation is actually constant and given by σ S . There is no dependency on the taxel distance D. For contact points close to the taxel and outside the taxel pair, the uncertainty of the force magnitude grows with distance (p 3 , p 4 ). The standard deviation for position localization σ P is typically largest in the center between the taxels and gets smaller on either side (exceptions are discussed in the following). The exact equations to compute these standard deviations are given in Sec. A and visualized Fig. 2 B, as described below.
Knowing the position accuracy σ P , we can quantify the super-resolution capabilities. For simplicity, we use 2σ P as the size of a virtual taxel, as shown in Fig. 2 A-I, corresponding to a confidence interval of about 68%. Thus, the spatial resolution is 2σ P and between two sensors at a distance D we can distinguish D/2σ P virtual taxels, which we define as super-resolution factor w.r.t. real taxel number n (n = 2 in Fig. 2 A-I): Influence of the isoline shape The shape of the TVIs depends on the properties of the transmission medium and the sensor unit type. We study the impact of the TVI shape on the accuracy of single contact force inference. We follow the general model that the response of a taxel to a force on the surface decreases monotonously with distance from the taxel center. This in turn leads to a monotonously increasing TVI.
We assume the attenuation of taxel response behaves as s ∝ 1/|d| α for a force at distance |d| from the taxel center with the attenuation exponent α. Without loss of generality, we consider The place where the lines corresponding to the particular sensor readings cross is the contact location that is marked as p 1 . The measurement noise σ S leads to uncertainty of position σ P and force σ F . A-II: Different intersection types. B: Effect of TVIs shape on super-resolution characteristics. B-I: TVIs for two taxels (at distance 1) with different attenuation exponents α.
Note the different shapes of the intersection areas. B-II: Resulting sensitivity, standard deviation of position localization and standard deviation of force inference. C: Theoretical super-resolution characteristics of a 1D sensor with multiple taxels. C-I: Three taxels localizing a single contact. C-II: The spatial distribution of accuracy for a single contact. Below the orange dash-dotted line, no super-resolution localization is possible. Notice the increase in accuracy for higher forces because multiple taxels are activated.
two taxels at a distance of 1. Their TVIs are given by where g(S) is the force corresponding to the measurement at distance 0. The TVIs for different attenuation behaviors are shown in Fig. 2 B. The accuracy of super-resolution sensing is strongly affected by α. For instance, for linear attenuation (α = 1), the localization will only be possible between the taxels, but not outside (unbounded intersection area). For concave curves (α < 1), there can be three intersection areas: one is between the two taxels and the other two are outside, which makes the reconstruction not unique. For the sake of comparison, we assume that in this case the correct intersection is known. Smaller attenuation exponent α yields a smaller sensitivity (bigger F S ) between the taxels but in general a more homogeneous distribution ( Fig. 2 B-II). The position accuracy, measured by the standard deviation, σ P is constant for α = 2. For normalized TVIs, as used here for comparison, σ P = σ S . The general case is given below. For α > 2 and α < 1 we find a reduced position accuracy between the taxels. Interestingly, for α > 2, the reduced position error at the outside is paid by a reduced force sensitivity. The ability to infer the force magnitude, given by σ F , is constant between the two taxels, irrespective of α. For contact points outside, larger values of α are better, whereas for α ≤ 1 no detection is possible. For 1 < α < 2, the position error and force error are both high outside. From this analysis, α = 2 yields the best overall characteristics, namely having high force sensitivity (small F S ), high position, and force accuracy (small σ P , σ F respectively). Values α > 2 trade off some super-resolution by sensitivity, but are generally also good. Taxel response attenuation with exponent α < 2 should be avoided, if possible.
An approximate relationship between σ P and σ S for the more general TVIs I S (d) = g(S) + λ|d| α is given by: where m 1 = λα|d| α−1 and m 2 = λα|D − d| α−1 are the absolute derivatives of the TVIs as detailed in Sec. A. Using the position accuracy at d = D /2 (which is the worst for typical exponents, α ≥ 2) for n = 2 sensors we can compute the super-resolution factor analytically using Eq. 4: Multiple taxels in a line A real sensor should clearly contain more than 2 taxels, so let us consider multiple taxels (with quadratic isolines) equidistantly placed in the transmission medium (as shown in Fig. 1 A-I). As expected, in the area where more than two taxels respond to the contact stimulation, higher localization accuracy is possible, as shown in Fig. 2 C-I. This leads to a reduction in uncertainty about the contact force. The intersection area is smaller and, due to averaging independent noise measurements, the variance also reduces with an increased number of TVIs intersecting. The sensitivity F S is shown in Fig. 2 C-II as a dash-dotted lower bound, together with the accuracy of localization σ P and force quantification σ F . The sensitivity is not homogeneous and is higher between taxels, which might be a surprising result at first glance. In summary, the most important take-home message is to have multiple taxels responding to a contact force because it improves accuracy. Closer placement of taxels increases both sensitivity and accuracy, whereas reducing taxel measurement noise improves accuracy but not sensitivity.
Multiple contact points In many applications, we are interested in detecting multiple simultaneous contact points. To detect two contact points, we need at first glance two pairs of taxels. However, when the contact points are too close, spurious contact points would be detected due to additional intersections of TVIs. The basic condition for successfully distinguishing them is shown in Fig. S2. There should be at least two taxels between these two contacts, and at most one taxel of them that is evoked by one contact not the other. With higher forces applied by these two contacts, the distinguishable distance for double simultaneous contacts is larger.

Physical Implementation in 1D
External contact at the sensing surface causes deformation of the transmission medium (TM) that can be measured by physical sensor units (taxels). On the one hand, we investigate the influences of taxel placement and TM (e.g., elastomer) properties on the TVIs and on the sensitivity. On the other hand, we augment our theoretical analysis by measuring the response curves of three suitable real physical sensors and choose one of them to validate our proposed theory.
Physical factors of TM influencing TVIs Where to place the sensor units within the transmission medium? How thick should be the transmission medium? What is the effect of material properties such as the Poisson's ratio? A static mechanical model simulated with the finite element method (FEM) using Ansys (29) is built to answer the above questions. Details are explained in Sec. C. It indicates that the dependency of the deformation on the depth has an impact on the TVIs, and suggests placing sensors more close to the sensing surface for higher sensitivity. Considering the structure design, thicker material causes more displacement and in turn increases sensitivity, because the fixed boundary effect becomes less dominant. From the material property perspective, the Young's modulus and the Poisson's ratio are the two main properties of the TM we consider. The Young's modulus describes how easy the material deforms and has a proportional impact on the deformation. A soft material (small Young's modulus) improves sensitivity but is also increasingly deformed by inertial effects, which needs a trade-off with the material density. The Poisson's ratio measures the relative transverse/axial expansion when the material is axially compressed. Most elastomers have Poisson's ratio around 0.5, and metals have around 0.3 (30). Decreasing the ratio, the radial displacement becomes much less sensitive, whereas the depth displacement has higher sensitivity (lower TVIs  TVIs of real sensors We consider strain gauges that measure the change in curvature along one direction averaged over the sensing area (16,17); accelerometers that are able to measure the absolute inclination of the local elastomer patch using the gravity direction as a reference (31); and barometers that sense the volume change caused by the material's deformation in the form of isotropic pressure (32, 33). As detailed in Sec. D, the strain gauge has a non-monotonic behavior where one strain gauge value has several position-force possibilities. The accelerometer does not have this problem but has a tiny "blind spot" directly above the sensor unit (no inclination), however super-resolution localization is well possible. The barometer shows the convex and monotonic properties as described in our proposed theory.
Quantitative analysis in 1D Due to the convexity and symmetry of the barometer's isolines, we use it to validate our proposed theory. We mold six barometers in the elastomer along a straight line with approximately 6.5 mm distance to each other as shown in  Fig. 3B. They suggest a very good super-resolution potential. As mentioned above, the measurement noise (σ S ) in the sensor units introduces uncertainties in the position and force strength inference (σ P and σ F ). To solve the inverse problem of predicting the indentation position and force magnitude from the sensor measurement, we employ a machine learning model using squared error loss, which yields a prediction with minimal variance. In this way, we circumvent a manual computation of intersection areas, which would pose problems with real-world deviations from the idealized TVIs. Nevertheless, our theory can be seen as an upper bound for the performance of the machine-learning performance. For the architecture and training details, see Sec. E. The results are summarized in Fig. 3 The inference accuracy of the position and force magnitude is higher with stronger indentation force, but generally very accurate (evaluated at locations that were not included during training), as shown in Fig. 3 C-I. Figure 3 C-II shows the averaged super-resolution factors Ω (Eq. 4) predicted by our theoretical analysis and achieved by the machine learning solution, which have a remarkable correspondence. As expected, the theoretical predictions are higher, as they represent the best achievable result. Higher force allows for higher super-resolution. Averaging it over the applied force range (from 2 mN to 1.5 N), we obtain average super-resolution factors of 109 with the ML method compared to the prediction of 165 from the theory. The spatial resolution of position error and force magnitude error of the prediction models are shown in Fig. 3D. To show the impact of different TVI overlaps, we present in Fig. 3 E-I a sensor with varying distance between the sensor units. Figure 3 E-II shows the spatial resolution of position error thereof, which reveals the direct effect of the distance between adjunct sensor units. In both cases, the overall shape resembles our theoretical prediction, cf., Fig. 2: higher errors occur in locations where less TVIs overlap.

Super-resolution in 2D
The above analysis for the 1D case helps us to investigate a sensor with a flat or curved 2D sensing surface. To simplify the analysis, we continue to assume a homogeneous transmission medium and an isotropic sensor unit, which is, for instance, approximately true for barometric sensors but violated for strain-gauge sensors. We consider a flat 2D sensing surface with coordinates x and y. The concept of isolines translates into isosurfaces, as shown in Fig. 4 A-I. However, we still call them TVIs for consistency. Clearly, with only two taxels, an accurate localization can only be done along one dimension (x) but not along two (x and y), as shown in Fig. 4 A-II. To make proper super-resolution localization possible, at least three taxels need to respond to a stimulus, as shown in Fig. 4 A-III. As before, the accuracy increases if more taxels are involved, and the sensitivity is higher between taxels over the sensing surface ( Fig. 4 A-III). As detailed in Sec. F, to localize two simultaneous contacts, at least 6 taxels are required to respond. Similar to the 1D case, if the contacts are too close, spurious intersections can occur. The spurious intersections can be ruled out because of high elongation in one direction, similar to Fig. 4 A-II. This is a new feature that was not observed in 1D. Very similar considerations are also valid for curved sensing surfaces, which is illustrated in Fig. S6 C-IV.
We mold 25 barometers in the elastomer with a 5 × 5 grid layout with approximately 6.5 mm distance to each other as shown in Fig. 4  In coherence with our theory, at the boundaries this improvement is smaller because less TVIs overlap. In comparison, the spatial distribution of the force accuracy is relatively homogeneous. The accuracy decreases slightly with increased indentation depth, and the boundaries tend to have an underestimation of force magnitudes. Notice, our testbed has a position precision of 0.05 mm, the force-torque sensor has a force precision of 0.01/0.01/0.02 N (F x /F y /F z ).

Discussion
We present a new way to characterize, analyze, and predict force sensation at super-resolution for haptic sensors. Our theory is based on sensor isolines that allow for a direct assessment of the uniqueness of contact position reconstruction. We derive quantities such as minimum force sensitivity, localization, and force sensing accuracy. These allow us to analytically compute the super-resolution factor, namely, the number of distinguishable locations between two real sensor units. With the help of an FEM model, we give guidelines for common design choices, such as placement of the sensor units within an elastomer as well as material properties. We analyze three commonly used sensor types: strain gauges, accelerometers, and barometers. Both, accelerometers and barometers, show the necessary properties for super-resolution sensation. We conduct two cases study to evaluate our theory using a line (1D) and a grid (2D) of barometer sensor units embedded in elastomer skin. Following our theoretically derived design choices enables us to achieve a remarkable performance of 109-fold/1 254-fold super-resolution for the 1D and the 2D layout using machine learning methods. Our theory predicts a super-resolution factor of up to 165 for the 1D case. We hope that our approach can help the design of new and capable haptic sensors. The major insights from our study are: The sensor units in the transmission medium (elastomer) should have convex isolines. It is beneficial to have the sensor units "float" in the center of the elastomer or closer to the sensing surface. Flexible wiring helps to have the best sensitivity. A thicker elastomer layer seems beneficial, and materials with small Young's modulus, high Poisson's ratio, and big yield strength are recommended. The distance between sensor units should be such that for a majority of forces, the isolines of neighboring taxels intersect. For inferring simultaneous contacts, a single sensor unit should be only activated by a single contact, and the distinguishable distance increases when the force magnitudes of multiple contacts are higher.
The actual information inference can be implemented using machine learning, as we present here. The decisive quantity for super-resolution perception is spatial distribution of the uncertainty of contact location and force magnitude inference. Exploiting the equivalence of minimizing model uncertainty by maximum likelihood estimation and least-squares error minimization for Gaussian residuals, the machine learning model directly optimizes for the quantities of interest. Whereas the theory describes an upper limit for the super-resolution capabilities under the model-assumptions, we find that the ML methods approach the predicted accuracy in the considered real-world sensors. In this way, our theory offers guidelines towards the system-level design of machine-learning-driven haptic sensors. Besides the above-mentioned mechano-electrical properties, the predictions by the theory can be used to validate the suitability of the employed ML model, the data collection, and the training procedure, as we can expect the model and the predictions to be aligned. For instance, inaccuracies in the testbed or recording a lot of data in regions without TVI overlaps can lead to suboptimal overall performance.
Future work will be devoted to analyze the accelerometers in more detail and investigate shear forces. Furthermore, structured transmission media e.g. with ridges (34) and with multiple layers (35) are an interesting direction.

Methods
We conducted several experiments to make informed design choices and validate the proposed theory.
Finite element analysis To analyze the physical factors of the transmission medium influencing TVIs, we built a suitable finite-element model using Ansys (29). Testbed We create a custom testbed with three degrees of freedom (DoF). Three DoF control the Cartesian movement of the probe ( x, y, z) using linear guide rails (Barch Motion) with a precision of 0.05 mm. The 4 mm spherical probe is fabricated from an aluminum alloy and is rigidly attached to the Cartesian gantry via an ATI Mini40 force/torque sensor with a force precision of 0.01/0.01/0.02 N (F x /F y /F z ). The sensor is positioned on the base, and the indenter is used to contact it at the desired location.
Data To obtain a variety of normal forces, the indenter is moved to a specified location, touches the sensing surface, deforms it increasingly by moving normal to the surface with fixed steps of 0.2 mm. After a pause of 2 seconds to allow transients to dissipate, we simultaneously record the contact location, the indenter contact force from the testbed's force sensor, and the barometric units' values from the molded sensors. All the data are collected and combined using a standard laptop. More specifically, for the 1D sensor (Fig. 3), the testbed makes the indenter contact 2501 positions evenly spread along the sensor center line (from -25 mm to 25 mm) with 40 incremental indentation depths (0.1 mm each) at each position. For the 2D sensor (Fig. 4), 69 × 69 positions (a grid with 0.5 mm apart from each other) evenly distributed on the sensing surface (34 mm × 34 mm) are probed by the indenter with 20 incremental indentation depths(0.2 mm each).
The data of the "Minerva" pattern (Diameter: 26 mm) is collected similarly for validating the super-resolution functionality. The pattern has a position resolution of 0.007 mm. The indenter probes the sensor at each position with 10 incremental steps (0.4 mm each with a transients pause of 2 seconds).
TVIs and related parameters We implement the following steps to compute the isolines ( Fig. 3 and Fig. S4). First, we linearly interpolate the sensor values and force values. Second, we choose a position-related sensor value and find the corresponding force measurement in that position. Third, we draw the position-force curve for that sensor value with the same color as shown in Fig. 3B & E. Based on the derived isolines, we use curve function (λ · d α + g(S)) to fit each derived isoline. The constant c is calculated by dividing g(S) over S. The super-resolution factor Ω is calculated thereby using Eq. 8.
Machine learning For the 1D sensor, we use a standard MLP (multi-layer-perceptron) with six fully connected hidden layers with 100 neurons each. The data consist of 50 k samples that are split into datasets of training, validation, and test with a ratio of 3:1:1. The data within the region of from -16.25 mm to 16.25 mm are selected as training data from the training dataset, and this region spans exactly from the center of the left barometric unit to the right one. The models are trained with the squared error loss, the Adam optimizer (learning rate: 5 · 10 −4 , epsilon: 10 −5 ), and a batch size of 200 samples in 1 million iterations. The models for position and force inferences are separately trained using the same architecture and training settings. For the 2D sensor, we use an MLP with ten fully connected hidden layers with 100 neurons each. The data consists of 95 k samples that are split into datasets of training, validation, and test with a ratio of 3:1:1. The models are trained with the L2 loss, the Adam optimizer (learning rate: 2 · 10 −4 , epsilon: 10 −5 ), and a batch size of 100 samples in 1 million iterations. The models for position and force inferences are separately trained using the same architecture and training settings.

Theory and Design of Super-resolution Haptic Skins
The PDF file includes: Page ii: Sec. A: Computation of σ P and σ F .
Page ii: Fig. S1 Computation of σ P and σ F .
Page iv: Sec. B: Multiple simultaneous contacts discrimination.
Page iii: Fig. S2 Multiple contacts discrimination.
Page v:

A Computation of σ P and σ F
We first go through the accurate computation of the inference uncertainties for the idealized TVIs that is used in the numeric results. Then we present an approximation that we use to compute the super-resolution factor analytically.
Precise computation of σ P and σ F For two taxels at distance D, as shown in Fig. S1 A, and the TVIs described by the following equations: we use the exact corners of the intersection area of the two TVIs with their measurement uncertainty hose. These intersection points are computed numerically, and we use this method for all plots in the main paper. In more detail, for each TVI we have an upper and lower bound for the uncertainty (one standard deviation) denoted by h ± : There are 4 intersection points: Fig. S1 B. The size of the bounding box of these four points defines σ P and σ S .  . When these two contacts happening simultaneously, the value of the sensing unit S-3 turns to be an orange-dash-dotted line (combined TVI). This will result in an unreliable inference. Multiple spurious contact points occur (" ") rather than two actual contact points ("•") Approximate analytical computation of σ P Here, we give a closed-form expression relating the force measurement noise to the position accuracy between two taxels. We make a first-order approximation of the intersection area as a parallelogram, as shown in Fig. S1 C. For two taxels at distance D and the TVIs described by Eq. S1 and Eq. S2, we can compute the (absolute) derivatives of the TVIs at the intersection point d: The trigonometric construction is shown in Fig. S1 C and yields the following result:

B Multiple Simultaneous Contacts Discrimination
In many applications, we are interested in detecting multiple simultaneous contact points. In this section, we investigate the criterion to discriminate double contacts. At first glance, we need 2 pairs of taxels. However, when the contact points are too close, spurious contact points would be detected due to additional intersections of TVIs. The basic condition for successfully distinguishing them is shown in Fig. S2. We consider two contact points with individual force strength are located in a 1D sensor layout with six sensing units. When one contact is located in the left green-shaded area, then a second contact can be clearly discriminated when it is within the right green-shaded area. In contrast, when the second contact is located in the gray-shaded area, it activates at least three sensing units (S-3,S-4,S-5). Notice, the sensing unit S-3 is also activated by the first contact. Thus, the value of S-3 is then not reliable for triangulation inference and will create spurious inference points and modify the position of the otherwise correctly inferred locations. The most important take-home message is to have at least two taxels between these two contacts and at most one taxel of them that is evoked by one contact not the other. With higher forces applied by these two contacts, the distinguishable distance for double simultaneous contact is larger.

C Physical Factors Influencing TVIs
We investigate the detailed influences of the properties of the transition medium (e.g. elastomer) and the taxel placement on the TVIs and on the sensitivity. External contact at the sensing surface causes deformation of the transmission medium that can be measured by physical sensor units (taxels). In this section, we assume the taxels are able to measure deformation/strain within the medium, either isotropically or directionally. In principle, the deformation can be described by the absolute movement of elements, called displacement, and by the relative movement of elements, denoted as strain. To study the isolines we consider the displacement as a measure of deformation. The displacement is computed for a static mechanical model simulated with the finite element method (FEM) using Ansys (29). The questions we want to answer in this section are: Where to place the sensor units within the transition medium? How thick should be the transition medium? What is the effect of material properties such as the Poisson's ratio?
Simulation model The full model is a cylinder-shaped sensor transmission medium with a diameter of 200 mm and a thickness of 5 mm, 10 mm, or 15 mm, but we simplify the model in the FEM simulation into a 2D plane due to the axis-symmetric property of the cylinder and the spherical indenter (10 mm diameter) as shown in Fig. S3 A-I. The transmission medium is constrained at the bottom by a fixed boundary. We analyze the displacement when a normal force is applied to the contact surface via the indenter. In Fig. S3 A-II the displacement is evaluated in xand z-direction, as well as in terms of total displacement. Considering the displacement v at three different depths, close to the sensing surface (L1), in the middle (L2), and close to the fixed boundary (L3), shows significant differences. The displacement in depth direction (∆z) decreases monotonously with distance, but is much stronger at L1 than at L3. The radial displacement (∆x) first increases and then attenuates, which will give rise to non-monotonous TVIs.

Taxel depth
The dependency of the deformation on the depth has an impact on the TVIs. Figure S3 B-I shows the TVIs for ideal displacement sensors at three different depths (L1-L3) described above. We plot the TVI for the smallest non-zero taxel value. However, because there can also be negative displacements, we add a light-shaded TVI for small negative values if occurring. When measuring depth deformation, the TVIs have a convex shape for positive sensor values and a second parabola at either side for the negative sensor values. This holds for all depths, although only visible for L1 due to a limited plot range. The radial direction is more surprising, as the center is not moving radially, the TVIs have a pole in the center. The total deformation is monotonous for low to middle taxel depths (L1 and L2). Depending on the physical quantity measured by the real sensor unit, this leads to different conclusions: For strain gauges and accelerometers, for instance, the global displacement can be considered as a good proxy, as they measure the curvature and inclination of the material. Thus, these sensors should be more close to the sensing surface. Barometric sensors, for instance, measure local displacement, for which we cannot give a recommendation at this point.
Sensing surface thickness Guided by the previous result, we consider a taxel at depth L1 (1 mm below the sensing surface) for different total thicknesses of the elastomer. Figure S3 B-II presents the TVIs. The radial deformation is, as expected, only a little effected by the amount of material underneath the taxel. The depth deformation, on the other hand, shows a much larger effect of thickness, because for smaller thickness we get less displacement. In summary, the thicker material causes more displacement and in turn increases the sensitivity, because the fixed boundary becomes less dominant. Notice, for thicker material, the shear/x-directional displacement has stronger attenuation property and has less shear sensitivity. Nevertheless, there are several ways to increase the sensitivity in the shear direction, as proposed in earlier work (34,35). One is to add ridges on the contact surface (34). Another is to use multiple layers of materials to cause mechanical contrast and affects the force distribution (35).
Material properties For a robust sensor design, we need to make sure that the contact-caused stress is smaller than the yield strength of the transmission medium, and the material deformation is an elastic behavior. The main two properties of the transmission material we consider are the Young's modulus and the Poisson's ratio. The Young's modulus describes how easy the material deforms and has a proportional impact on the deformation. A soft material (small Young's modulus) improves sensitivity, but is also increasingly deformed by inertial effects. Poisson's ratio measures the relative transversal/radial expansion when the material is axially vi compressed. Most elastomers have Poisson's ratio around 0.5 and metals have around 0.3 (30). Figure S3 B-III shows the TVIs for different Poisson's ratios. Decreasing the ratio, the radial displacement becomes much less sensitive, whereas the depth displacement has higher sensitivity (lower TVIs). Thus, depending on the measurement direction of the real sensor unit, different Poisson's ratios are preferred.
Indenter shape What is the impact of the indenter shape (or the object getting in touch with the sensing device later) on the TVIs? For illustration, we simulate the displacement caused by three different indenter shapes: a sphere, an ellipsoid, and a cylinder. The displacement decreases and becomes flatter for the cylinder near the contact center, as shown in Fig. S3 B-IV. The cylinder with flatter contact surface causes less sensitivity for shear displacement. Increasing the indenter size, the displacement becomes smaller, which results in less sensitivity too. As a remark, for too small indenters, there is a risk to exceed the material yield strength and break the contact surface material.

D Taxel Value Isolines of Real Sensors
We complement our theoretical analysis by measuring the response curves of three suitable sensors. We consider strain gauges that measure the change in curvature along one direction averaged over the sensing area (16,17); accelerometers that are able to measure the absolute inclination of the local elastomer patch using the gravity direction as a reference (31); and barometers that sense the volume change caused by the material's deformation in the form of isotropic pressure (32, 33). Fig. S4 A, we mold a single sensor unit at the center of an elastomer disc (120 mm diameter and 10 mm thickness) and measure the sensor unit's response when indenting it with a controlled depth along two perpendicular axes using an automated 3-DoF testbed with a force sensor. The strain gauge sensor is EP08-250AE-350, the accelerometer is ADXL345, and the barometer is MPL3115A2. We solder extra thin wires (CU-enameled wire with a diameter of 0.15 mm, ME-MeßSysteme GmbH) to the chips to be able to mold them inside the elastomer with minimal mechanical influence of the taxel. We mold these three sensors floating in the middle and center of three elastomer pieces (EcoFlex 00-30) individually. Additionally, we build a 3-DoF test bed (Barch Motion, precision: 0.05 mm) equipped with a 3-axis force measurement device (ATI Mini40, precision: 0.01/0.01/0.02 N (F x /F y /F z )) to indent the elastomer surface in a precise and automated way. Data processing Based on the recorded data, we present in Fig. S4 B the sensor values, the applied forces, and the TVIs for each sensor type as a spatially changing quantity. We only show the data for one of the two directions. For the strain gauge unit, this is along its measurement direction. Accelerometers are interesting, because they can distinguish between both directions. The barometer is fairly isotropic up to small deviations because of the rectangular sensor housing. We implement the following steps to compute the isolines: first, we linearly interpolate the sensor values and force values; second, we choose a position-related sensor value and find the corresponding force measurement in that position; third, we draw the position-force curve for that sensor value with the same color. The strain gauge has a non-monotonic behavior, where one strain gauge value has several position-force possibilities. The accelerometer does not have this problem, however, it has a "blind spot" directly above the sensor unit (no inclination), however, this is a tiny area. Note that the TVIs on both sides are for different sensor values, so super-resolution localization is well possible. We leave a theoretical analysis of these TVI shapes for future work. The barometer shows the convex and monotonic properties as described in our proposed theory. Notice the two orange "wings" which are the TVI for negative values caused by the lateral pulling force from distant contacts.

Experimental setup As shown in
TVI fitting and sensor-force relationship (constant c) To support the assumption for our theory introduced in Section "The Model", we herein introduce the analytical steps for approximating taxel-value-isolines (TVIs) from the sensor and force measurements. We first empirically acquire the sensor isolines, and then fit a function to each individual isoline, as shown in Fig. S4 C-I. Taking the g(S) of each barometric values at d = 0 as the force strength, we can acquire the relationship between sensor value and force strength captured by the constant c (see text below Eq. 2). The value of c is simply the ratio between the applied force g(S) and the barometer value. As shown in Fig. S4C-II, the linearity assumption holds. Moreover, the attenuation coefficient λ and power α are also approximated. As shown in Fig. S4C-III, λ becomes larger with increased barometer value, while α attenuates approaching α = 2.

E 1D Sensor
Theory Due to the convexity and symmetry of the barometer's isolines, we use it to validate our proposed theory. We mold six barometers in the elastomer along a straight line with approximately 6.5 mm distance to each other as shown in Fig. 3 A. As mentioned above, the measurement noise (σ S ) in the sensor units introduces uncertainties in the position and force x strength inference (σ P and σ F ). There are two major noise sources in our setting, one is the noise measurement in the six barometers σ Si , another is the force-torque sensor σ SF . We collect 1000 samples of the force-torque sensor and the barometer values in the unloaded case to evaluate the measurement noise levels, as shown in Fig. S5A  : to calculate the super-resolution factor Ω i for each barometric sensor unit, assuming two identical units are used. The overall Ω of the whole 1D sensor is 165 that is averaged over these six Ω i for each discrete threshold of 0.02 N over a range of (0.02 N-1.5 N).
Machine Learning To solve the inverse problem of predicting the indentation position and force magnitude from the sensor measurement, we employ a machine learning model using squared error loss, which yields a prediction with minimal variance. In this way, we circumvent a manual computation of intersection areas, which would pose problems with real-world deviations from the idealized TVIs. Nevertheless, our theory can be seen as an upper bound for the performance of the machine-learning performance. We use a standard MLP (multi-layer-perceptron) with six fully connected hidden layers with 100 neurons each. The data consist of 50 k samples that are split into datasets of training, validation, and test with a ratio of 3:1:1. The data within the region of from -16.25 mm to 16.25 mm are selected as training data from the training dataset, and this region spans exactly from the center of the left barometric unit to the right one.

F 2D Sensor
Theory The above analysis for the 1-dimensional case helps us to investigate a sensor with a flat or curved 2-dimensional sensing surface. To simplify the analysis, we continue to assume a homogeneous transmission medium and an isotropic sensor unit, which is, for instance, approximately true for barometric sensors but violated for strain-gauge sensors. We first consider a flat 2D sensing surface with coordinates x and y. The concept of isolines translates into isosurfaces, as shown in Fig. S6 A-I. However, we still call them TVIs for consistency. Clearly, with only two taxels, an accurate localization can only be done along one dimension, but not along two. The intersection volume of the TVIs is presented in Fig. S6 A-II & III. For attenuation exponents, α = 2 the intersection volume changes depending on the transversal position y as shown in Fig. S6 B. To make proper super-resolution localization possible, at least 3 taxels need to respond to a stimulus, as shown in Fig. S6 C-I. As before, the accuracy increases if more taxels xi are involved, as shown in Fig. S6 C-II for 4 taxels arranged in a regular grid. We also compute the sensitivity distribution over the sensing surface for both sensor arrangements. The honeycomb pattern allows for a more homogeneous sensitivity. To detect two simultaneous contacts, at least 6 taxels are required to respond. Similar to the 1D case, if the contact points are too close, spurious intersections can occur. In some cases, as shown in Fig. S6 C-III, spurious intersections can be ruled out because of high elongation in one direction. This is a new feature that was not observed in 1D. Very similar considerations are also valid for curved sensing surfaces, which is illustrated in Fig. S6C-IV.
Hardware and Machine Learning We mold 25 barometers in the elastomer with a 5 × 5 grid layout with approximately 6.5 mm distance to each other as shown in Fig. 4 Fig. S7, which show approximately isotropic properties along all directions. We employ two machine learning models to predict the indentation position and force magnitude, respectively. We use an MLP with ten fully connected hidden layers with 100 neurons each. The data consists of 95 k samples that are split into datasets of training, validation, and test with a ratio of 3:1:1. The models for position and force inferences are separately trained using the same architecture and training settings.

Super-resolution
Factor Ω Using the trained machine-learning models, we evaluate the superresolution factor of the 2D sensor as the fraction between the number of virtual taxels n v and the real number of sensors (taxels) n r = 25 with the following equation: where A = 26 mm ·26 mm is the center sensing area of interest, and A v = π · σ Px · σ Py is the virtual taxel area calculated as an ellipse with radii of the position error (standard deviations) in x and y directions.  Fig. S6. Taxel values isolines for a 2D sensing surface. A-I: Taxel value isoline (iso-surface) for a single taxel. A-II & III: Intersection volume (due to measurement uncertainty) for two taxels at a distance D along the x axis, see top view. The localization would be very uncertain along the y direction. B: Intersection volume for different attenuation exponents α. C: Proper localization requires at least 3 taxels for 1 contact point. I shows a hexagonal sensor placement with resulting sensitivity over the surface and an illustration of the intersection volume of the TVIs from 3 taxels. II is the same for a grid and 4 taxels. III illustrates a spurious contact localization for 2 contact points and 6 taxels (marked with "?"), however its uncertainty would be very large. IV illustrates the intersection of TVIs for the case of a curved sensing surface. The same applies as in the case of planar 2D, except that distances need to measured as geodesics on the curved surface.