> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <1

Deformation Matching: Force Computation
based on Deformation Optimization*
Koichi Hirota, Yusuke Ujitoko, Sho Sakurai, Takuya Nojima, Member, IEEE
Abstract—The tactile information to be presented to a user
during interaction with a virtual object is calculated by simulating
the contact between the object model and user model. In thesimulation, a distributed force is applied to the contact area on theskin tissue of users’ hands and results in deformation of the skintissue. The skin deformation caused by the distributed force is thetarget contact state that should be presented by the device.
However, most multipoint haptic displays do not have sufficientdegrees of freedom (DoF) to represent the target contact state.
This paper presents the concept and formulation of “deformation
matching,” whereby the output force is calculated to minimize theerror between the target skin deformation and skin deformationthat can be realized by the limited DoF device’s output force. For
comparison, the conventional concept of “force matching” wasalso formulated. The difference in human perception between
these two concepts in the expression of friction was investigatedthrough experiments using a pin-array tactile display capable ofstimulating 128 points. It was demonstrated that the perception of
the friction coefficient was more sensitive and the perception of thefriction direction was more accurate in deformation matching
than in force matching.
Index Terms—Haptic interface, Haptic rendering, Tactile
display, Virtual reality.

I. INTRODUCTION
Increasing the degrees of freedom (DoF) and presenting
distributed stimuli is a development trend of haptic devices. It
is expected that increased DoF will make it possible to convey
the local shape and surface properties of the object. Tactile
sense is a measurement system that is capable of receiving
distributed stimuli on the skin surface, and haptic devices needto have sufficient spatial resolution and DoF for the properties
of such receptors. Fingertips are known to have particularlyhigh tactile resolution, and tactile displays that exert force on
multiple points, hereinafter referred to as multipoint tactile
displays, have been developed.
The tactile information to be presented to the user during
interaction with a virtual object is calculated by simulating thecontact between the object model and user model. Fingers andhands have flexible tissue on the surface, which exerts a
distributed force on the contact area with an object and resulting
in the deformation of the tissue. This corresponds to the target
Submitted for review on April 11, 2021.
This work was supported by JSPS KAKENHI (Grant Number JP17H00754).
In part, this research used the ADVENTURE System from the University ofTokyo. Moreover, this work was performed using body geometry data from the
BodyParts3D database provided by the National Bioscience Database Center.
K. Hirota is with the University of Electro-Communications, Chōfu, Tokyo
182-8585, Japan (e-mail: hirota@vogue.is.uec.ac.jp).

contact state that is presented by the device. The DoFs of thecalculated tactile information vary depending on the model
approximation degree. When the DoF of the simulation anddevice are equivalent, it is possible for the device to present thetarget contact state as it is. However, if the DoF differ, it isnecessary to calculate a suitable output adaptively for thelimited DoF device to realize the target contact state.
The aim of this study was to establish a haptic rendering
algorithm that calculates the output force of a multipoint tactile
display from the target contact state. The majority of studies on
multipoint tactile displays have focused on device
implementation, and the evaluation of such devices has been
performed by focusing on the benefits derived from the densityof the stimuli, such as improvement in the recognition of fine
two-dimensional patterns [1, 2, 3, 4].
We previously attempted to present the contact force byusing a pin-array tactile display, by approximately mapping theforce to closely located pins. However, in the preliminary
experiment, it appeared that it was difficult to provide a
sufficient sense of reality. One of the reasons for this is that therendering method that maps the force to the pin cannot
sufficiently reproduce the deformed state of the skin.
Recently, Perez et al. proposed the concept of ‘optimizationbased haptic rendering’ [5], which formulates the calculation of
the device output as an optimization problem that minimizes theerror between the target deformation and the deformation
expected to be caused by the device output. In this pioneering
work, an approach for minimizing skin deformation error has
been presented and demonstrated. In subsequent research, an
extension of the approach to optimization regarding skin stresshas been investigated [6]. However, their work only addressedthe 3 DoF device and did not address the higher DoF device.
In the current study, we focus on the deformation, rather thanforce, of the skin contact state. We propose a rendering method
based on the concept of minimizing the error between thedisplacement of the target deformation and that of the skincaused by the device. The device output force is calculated bysolving an optimization problem that considers the device
constraints. This approach is referred to as “deformation
matching” in this study.
The contributions of this study are the formulation ofY. Ujitoko is with NTT Communication Science Laboratories, Atsugi,
Kanagawa 243-0198, Japan (e-mail: yusuke.ujitoko@gmail.com).
S. Sakurai is with the University of Electro-Communications, Chōfu, Tokyo
182-8585, Japan (e-mail: sho@vogue.is.uec.ac.jp).
T. Nojima is with the University of Electro-Communications, Chōfu, Tokyo
182-8585, Japan (e-mail: tnojima@nojilab.org).

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <deformation matching and the evaluation of the effectiveness of
the approach through experiments using subjects. In theformulation, it is assumed that the relationship between theforce and skin displacement can be expressed by a deformation
model. Moreover, it is assumed that the device applies force,
rather than displacement, on the skin. Although the nonlinearity
associated with skin deformation has been noted in past studies[7], a linear model is used as the first step owing to the stabilityof convergence, as described below. Furthermore, as acomparative approach for the proposed method, we formulate
the idea of minimizing the error between the target force andacting force of the device, hereinafter referred to as “force
matching.” This method is considered to be a generalization of
the conventional method of calculating the device output by
distributing the target force to nearby pins. In the evaluation
using subjects, the effectiveness of the proposed method for thepresentation of frictional force is verified by means of a pinarray-type tactile display that exerts forces on the skin throughpins. The presentation of frictional forces is a common tactile
expression and is also important in applications; however, it isa challenging task for pin-array devices that can only exert force
in a direction that is approximately perpendicular to the skin.
We attempt to reproduce the deformation of the finger when
friction is applied by deformation matching.
The following section describes related research and clarifies
the focus of our study. The formulations of the deformation
matching and force matching are presented in Section III, andthe implementation thereof for the experiments is described in
Section IV. The conditions and procedures of the experiments
are detailed in Section V, and the results are discussed in
Section VI. Finally, Section VII summarizes the results of thestudy and outlines future issues.
II. RELATED RESEARCH
A. Haptic Devices
The need for haptic feedback was noted in the early stages ofresearch on tele-robotics and virtual reality (VR), and various
implementation methods have been proposed. These methods
can be categorized into three typical approaches. The first is anapproach that expresses the force acting on a tool. Many
commercial haptic devices (e.g., PHANToM with pen grips andFalcon) are based on this method. The second is an approachthat expresses the object itself or its properties. This method
provides the reality of an object by expressing values such as
the object shape, the forces associated with the operation (e.g.,
movement, inertia, and vibration), and surface properties (e.g.,friction and smoothness). The third is an approach that
expresses the tactile stimulus on the skin surface of the user.
This method creates a tactile experience by presenting theforces and displacements on the skin surface.
Various stimulus types have been investigated to express
tactile stimuli on the skin surface of the user. The fundamental
approach involves presenting the force. The original
PHANToM presented a force vector to a finger via a gimbal andthimble [8]. Moreover, the presentation of forces with two DoFs
(i.e., compression and shear) by using a wearable device has

2been investigated [9]. Methods that densely integrate actuators
have been developed to present a distributed force on the skin.
Furthermore, the use of gel materials [10], solenoids [11],flexible tendons [1], shape memory alloys (SMAs) [12], andpneumatics [2] has been proposed.
Another fundamental approach is the application of skindeformation. This method approximates the deformed shapewith a plane [13] and attempts to increase the DoF to improve
deformation accuracy [3, 4]. Skin stretch is another concept thatis similar to the approach of applying deformation to the skin.
It has been demonstrated that it is possible to present torque tothe palm [14], present action force to the fingers [15], andexpress softness [16]. The application of vibration stimulation
to the presentation of two-dimensional patterns such as
characters has been studied for a long time [17], and in recent
years, the application to texture expression in VR [18] has also
been studied. Electrical stimulation is a method of electricalaction upon mechanoreceptors and nerves. A method has beendeveloped for selectively stimulating the receptors bycontrolling an electrode array [19], which is expected to be
applied to general presentations such as pressure and vibration.
B. Haptic Rendering
In addition to the device, the calculation process of haptic
information is indispensable for haptic feedback. This process
is known as haptic rendering. According to Salisbury, “haptic
rendering is the process of computing and generating forces in
response to user interactions with virtual objects” [20]. Early
haptic rendering research investigated models for user–object
interaction through points of action, and in this context, Godobject methods [21] as well as shape and texture representation
methods [22, 23] have been proposed. Subsequently, a modelthat approximates the surface contact by the skin was examinedbased on the demand for the reality of the gripping operation by
the finger or hand. Moreover, an approximate deformation
model by a spherical surface [24], an approximate calculationof the contact area [25], a reaction force calculation by skinning[26], and aggregation of constraints by contact [27] have been
proposed. In recent years, techniques for physically computing
the skin deformation have been investigated, and theintroduction of locally deformable pads into finger models [28]and interactions with deformable models of the entire hand [29,30] have been studied. Furthermore, it has been noted that tissue
nonlinearity needs to be introduced for more accurate contact
state calculations [7].
Owing to increasing DoF and complexity of haptic devices,
the relationship between the target contact state and the contact
state presented by the device has become less trivial. Therefore,
a method of mapping from the former to the latter is required.
Sato et al. proposed and verified the concept of calculating theelectrical stimulation from the strain energy density associatedwith the tissue deformation resulting from contact for thepresentation of electrical tactile stimulation [31]. As describedlater, the firing frequency of the receptor is proportional to thestrain energy density. Maemori et al. also proposed a method tocontrol the suction display by using strain energy density as anindex [32].

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <In general, the contact state on the skin is not necessarily
directly controllable depending on the configuration of thedevice. In such cases, the device output that will realize or
approximate the target state must be estimated based on a modelthat relates the device output with the contact state. Hence, theestimation can be formulated as a problem for finding thedevice output that minimizes the error in the contact state using
the model. Such approaches are categorized as "optimizationbased haptic rendering".
A pioneering work has been done by Perez et al. in which
they proposed a method for calculating device configuration in
the presentation of a plane to a finger using a thimble tactile
device [5]. This technique was formulated as an optimization
problem that minimized the difference between the target skindisplacement and the displacement caused by the device on theskin with postural constraints of the device. Verschoor et al.
proposed and evaluated a method for matching the stressgenerated by the device on the skin with the target stress for a
thimble tactile device [6]. The stress caused by the device wasestimated using a skin deformation model. In this method, thefriction state was estimated considering the dependence on thedevice operation trajectory, and a neural network was used forthe calculation to obtain the stress distribution from the device
configuration and friction state.
Another research applies the concept of optimization-based
haptic rendering to the presentation of distributed pressureusing a mid-air ultrasound device. This research proposed andevaluated algorithms for finding an optimal path of thestimulation under constraints due to the physical characteristicsof the device and temporal resolution of tactile perception [33,34].
C. Tactile Perception
A typical receptive property of mechanical tactile stimuli isthe tactile two-point threshold, which allows humans to
distinguish between two points on the skin when they are
stimulated simultaneously [35]. The tactile two-point thresholdis often described as an indicator of the spatial resolution that isrequired for tactile displays. In fact, it has been noted thatrecognition with a spatial resolution higher than the tactile twopoint threshold is possible [36, 37]. Regarding the perception oftactile stimulus movement, it has been demonstrated that theinter-stimulus onset interval and duration affect the occurrenceof tactile apparent motion [38]. This suggests that it is necessary
to increase the stimulation density in the expression of smoothtactile movement.
Furthermore, research is currently being conducted toelucidate the perception of tactile sensation from a mechanicalperspective. Tactile receptor response properties have been
revealed to be consistent with strain energy density-basedmodels [39]. Moreover, it has been demonstrated that theplacement of receptors on the skin is consistent with theproperties of the receptors [40]. The illusion that skin suction
produces a sensation that is similar to skin compression has also
been reported. This illusion can be explained by the similarityof the strain energy densities between suction and compression
in the tissues near the receptors [41].

3

D. Focus of this Study
According to the above review, a high-density and multi-DoF
device is required to express the stimulus acting on the skinaccurately. However, a haptic rendering method for such a
device that considers the DoF of the actuator and the restrictionof the operation has not yet been established.
In this study, we propose a rendering method for a tactile
device that presents the force on multiple points on the skin.
The fundamental concept is the optimization of the device
output force based on the error of the skin deformation by thedevice from the target deformation. The minimization of theskin deformation error was also proposed in the study by Perezet al. [5], but their study focused on thimble tactile devices, andthe formulation of a multipoint tactile display was not discussed.
In this study, we prioritized computational stability, because
this work is a preliminary attempt at rendering multipoint tactile
displays. As an optimization index, optimization based on thestrain energy density near the sensory receptor, instead of thedeformation, was an option, but it was not adopted in this study.
This is because the illusion of suction [41] suggests that thisoptimization problem may have multiple local solutions. The
nonlinearity of skin deformation was also not introduced owing
to concerns that it would degrade the convergence of the
optimization calculations. Furthermore, this study employed a
pin array device that was previously developed in the laboratoryfor evaluation, and in this sense, the results described below is
specific to the device.
III. PROBLEM FORMULATIONIn this section, the deformation matching and force matchingapproaches are formulated. Deformation matching and force
matching are two different methods for calculating the device
output to minimize displacement and force errors, respectively,from the target deformation state.
A. Deformation Matching
Suppose that the relationship between the skin force 𝑓 anddisplacement 𝑢 can be calculated by the skin deformation
model as 𝑓 𝐾 𝑢 . Furthermore, it is assumed that the force
generated on the skin by the device action 𝑓 can be calculatedusing a device model as 𝑓 ∗ 𝐴 𝑓 . The deformation 𝑢∗ due
to the device action is expressed as follows:
𝐾 𝐴 𝑓 .
1𝑢∗ 𝐾 𝑓 ∗Deformation matching is the problem of determining 𝑓 tosatisfy the following:‖𝑢 𝑢∗ ‖ → 𝑚𝑖𝑛.
2
In the following, the formulation of a linear model isdiscussed. A model of skin deformation by the linear FEM isgenerally defined as follows [42]: The subscript c indicates the
DoF with fixed (or imposed displacement) boundary conditions.
For a given 𝑢 and 𝑓, 𝑢 and 𝑓 are calculated.
𝑢
𝑓
𝐾 𝐾
3
𝑢 .
𝐾 𝐾
𝑓
The device is operated by actuators (e.g., voice coil motors)
and presents an output force on the skin. The driving force of
the actuators is converted into an output force through a

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <transmitting mechanism that depends on the device. In thefollowing, we consider the case in which a linear relationshipexists between the driving force and the force acting on the skin.
Assuming that a device with 𝑛 DoF operates with a driving
force 𝑓 ( 𝑛1 ) to generate a force 𝑓 ∗ on the skin, therelationship between the two can be described as follows:4𝑓 ∗ 𝐴𝑓 .
In the pin-array device used in the experiment described later,
the acting force acts directly on the skin, without requiring atransmission mechanism, and the direction is limited to the pin
driving direction. Moreover, because the point of action of thepin and the model node do not always match, it is necessary to
associate the acting force of the pin with the neighboring nodes
according to the position of the point of action. Considering thefeatures of a pin-array tactile display, 𝐴 is defined as the
product of 𝐴 and 𝐴 (𝐴 𝐴 𝐴 ), where 𝐴 (3𝑛𝑛 ) is amatrix in which the vectors of the pin action directions arearranged, and 𝐴 (3𝑛 3𝑛 ) is a matrix of coefficients that
distribute the force of the pins to the model nodes.
In a linear model, the displacement 𝑢∗ due to the device
action 𝑓 is calculated by
5𝑢∗ 𝐾 𝐴𝑓 .
Let us consider the determination of 𝑓 that minimizes thedeformation error ‖𝑢 𝑢∗ ‖ for a given target displacement 𝑢.
‖𝑢 𝑢∗ ‖
𝑢 𝐾 𝐴𝑓→ 𝑚𝑖𝑛.
6
Note that the model nodes are not necessarily arranged in a
spatially even density. That is, the error defined by the equation
can be considered as an index that is weighted by the spatial
density of the nodes.
In the following, we define 𝐿 as 𝐾 𝐴 ≡ 𝐿 (3𝑛 𝑛 ).
‖𝑢 𝑢∗ ‖
𝑢 𝐿𝑓
𝑢 𝐿𝑓
𝐿𝑓
𝐿𝑓
𝐿𝑓 𝑢 𝑢 𝐿𝑓
𝑢 𝑢
𝑓 𝐿 𝐿𝑓 .
7
𝑢 𝑢 2𝑢 𝐿𝑓
Because the target displacement 𝑢 is a fixed value, 𝑢 𝑢 is aconstant term, and thus, the following minimization problem
arises:2𝑢 𝐿𝑓 → 𝑚𝑖𝑛.
8
𝑓 𝐿 𝐿𝑓In this case, because 𝑥 𝐿 𝐿𝑥
𝐿𝑥 𝐿𝑥 0 , 𝐿 𝐿 is asemipositive definite matrix. In reality, in a device without a0, 𝐿 𝐿 becomes a definite
control state where 𝑢∗ 0 for 𝑓matrix, and the above problem becomes a convex quadratic
programming problem [43]. In the unconstrained case, thissolution is provided using 𝐿 , which is the pseudo-inverse
matrix of 𝐿. This solution is referred to as DM0.
𝑓
𝐿 𝑢.
9
In many multipoint tactile displays, the point of action (orpin) and skin do not adhere to one another, and a force cannot
be exerted in the skin pulling direction. The pin-array device
used in the experiments described below has similar restrictions
because air pressure cannot be controlled by negative pressure.
This constraint is expressed in the following form:
0 𝑖 1, ⋯ , 𝑛 .
10𝑓
A straightforward means of introducing this constraint is toclamp the above solution (DM0) to 0. This solution is referred
to as DM1.

4

𝑓 ← max 𝑓 , 0 .
11
It is necessary to consider optimization under constraints toobtain a more precise solution. The necessary and sufficientcondition for the optimality of the conditional convex quadratic
programming problem is known as the Karush–Kuhn–Tucker(KKT) condition. The KKT condition in the above problem isexpressed as follows by using the Lagrange variable 𝜆 (𝑛1).
This solution is referred to as DM2.
2𝐿 𝑢 𝜆
⎧2𝐿 𝐿𝑓⎪𝜆𝑓
0
𝑖 1, ⋯ , 𝑛 .
12
𝑓
0
⎨
⎪
𝜆
0
⎩
B. Force Matching
As an approach for comparing the proposed methods in theprevious subsection, we also formulate a conventionalpresentation method focusing on force. The majority of
research on multipoint tactile displays has focused on the device,with few discussions on the calculation algorithm of the output
force. The general concept in this research context is toreproduce the calculated force on the device. This concept canbe generalized as force-based optimization; that is, the problem
of determining 𝑓 that minimizes the error between the target
force 𝑓 and force 𝑓 ∗ owing to the device action.
‖𝑓 𝑓 ∗ ‖‖𝑓 𝐴𝑓 ‖ → 𝑚𝑖𝑛.
13
As with DM0, the error defined by ‖𝑓 𝑓 ∗ ‖ is an index
that is weighted by the density of the nodes. Moreover, in thisdefinition, all the nodes are subject to calculation. The sum of
the forces of the nodes that are not in contact with the object oron which no external force acts calculated by the FEM is 0.
Therefore, this definition only considers the acting force of thenodes that are in contact with the object.
The formula for the force error is expanded as follows:‖𝑓 𝑓 ∗ ‖𝑓 𝐴𝑓
𝑓 𝐴𝑓
𝑓 𝐴 𝐴𝑓 .
14
𝑓 𝑓 2𝑓 𝐴𝑓
In the above, because the target force 𝑓 is a fixed value, 𝑓 𝑓
is a constant term. Thus, the following minimization problem
can be obtained:
2𝑓 𝐴𝑓 → 𝑚𝑖𝑛,15
𝑓 𝐴 𝐴𝑓where 𝐴 𝐴 is a definite matrix; hence, minimization issolved as a convex quadratic programming problem. Thissolution, without constraints, is obtained using 𝐴 , which is thepseudo-inverse matrix of 𝐴. This solution is referred to as FM0.
𝐴 𝑓.
16
𝑓
An easy means of considering the constraints of the force
feedback device is to clamp the above solution (FM0) to 0. Thissolution is referred to as FM1.
17𝑓 ← max 𝑓 , 0 .
It is necessary to consider optimization under the constraints
to obtain a more precise solution. The KKT condition of thisproblem is expressed as follows, using the Lagrange variable 𝜆1). This solution is referred to as FM2.
(𝑛

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <2𝐴 𝑢⎧2𝐴 𝐴𝑓
⎪
0
𝜆𝑓
0
𝑓
⎨
⎪
0
𝜆
⎩

5

𝜆
𝑖

1, ⋯ , 𝑛 .

18

IV. EXPERIMENTAL SETUP
The experiment focused on the expression of friction when
the finger touched a plane and the plane moved relative to thefinger. The rendering calculations were performed in two steps.
The first step was a contact simulation, which calculated thedeformation of the finger and the force acting on the contact
area due to frictional contact. This was the target contact statein the subsequent matching process. In the second step, thedevice acting force was calculated by solving the optimization
problem by means of deformation matching and force matching,as discussed in Section III. The same finger model was used forboth the contact simulation and matching. A pin-array device
that presented a pressure sensation at 128 points was used as aforce feedback device. Both the contact simulation andmatching were performed offline to generate time-series data
for the device output. In the experiment, the force was presented
to the subject by driving the device based on these data.
A. Simulation Model
The finger model represented soft tissue in the area beyond
the center of the intermediate phalanx (Fig. 1). The basic model
shape was obtained by cutting out a bone and skin shaped modelof the finger part from a body model database (BodyParts3d[44]). A tetrahedral mesh was generated for this shape using the
TetMesh program (ADVENTURE project [45]). Fixedboundary conditions were set for the nodes on the surface that
contacted the bones and nail. The mesh consisted of 27,737
tetrahedrons and 5,800 nodes, and the number of free nodes (n)was 4,435. The material was uniform in the mesh, the Young's
modulus was 𝐸
1.0 10 Pa, and the Poisson’s ratio was𝜈 0.48. The Young's modulus was adjusted at the time ofpresentation by using the calibration method described later.
In the contact simulation, the frictional force was generated
by the relative movement of the plane in contact with the finger.
In the initial state (time 𝑡 0), the finger was pressed against
the plane at a certain depth 𝑑. The simulation was performed
under four different conditions with 𝑑 values of 0.5, 1.0, 1.5,
and 2.0 mm. Over time, the plane underwent tangential
displacement 𝑙, as defined by the following function:𝑙
𝑅 cos 𝜔𝑡 .
19
In the above, the constant of the stroke is 𝑅 0.02 m and theangular velocity is 𝜔 𝜋/4 rad/s. That is, the plane
reciprocates from a displacement of -0.02 m to 0.02 m with aperiod of 8 s. The initial contact state at 𝑡 = 0 was calculated byapplying normal displacement of the surface and then waitingfor 100 simulation steps or 0.42 s for the simulation to become
stable. The frictional force acting on the finger was stable andunchanging shortly after the entire contact area started to slip.
In the experiment described later, it was assumed that a force
was presented for 2 s from the start of the tangential motion.
B. Device Model

Figure 1. Mesh model for contact simulation and matching. The soft tissue
in the area beyond the center of the intermediate phalanx was represented
by a tetrahedral mesh model. A fixed boundary condition was implemented
on the surfaces of the bones and nail, which is represented by the red mesh.
The pin-array tactile display was assumed to be mounted on the finger at asmall pitch angle. The points of action and lines of motion of the pins are
represented by the green dots and yellow lines, respectively.

Figure 2. High-density pin-array haptic display. The device has 128 pins,
the output force of which is controlled independently.

The experiment used a pin-array tactile display (see [2] fordetails). The device had 128 contact pins, and their output force
was controlled independently. The pins were driven
pneumatically, and the pin pitch was approximately 1.5 mm.
The pins were arranged on a cylindrical surface, and the driving
direction was toward the axis of the cylindrical surface (Fig. 2).
The pin stroke is approximately 5 mm, the maximum force is
0.4 N, and the control delay is approximately 75 ms. The weightof the device is almost balanced by the rubber band suspension,and the device is fixed to the finger by pressing the back of thefinger with a velcro fastener. The point of action of each pin on
the mesh model was calculated as the intersection of the pin
motion axis and the model surface without deformation when
the device was attached to the finger pad. The direction of theforce action by the pin was assumed to be equal to the directionof the pin operating axis (green dots and yellow lines,
respectively, in Fig. 1). In practice, the finger and pin may slidein a direction orthogonal to the motion axis; however, in themodel in this study, it was assumed that the point of action did
not change. Similarly, the friction between the pin and finger
may generate a force in the direction orthogonal to the
movement; however, this was also not considered.
As mentioned in Section III.A, 𝐴 ( 3𝑛𝑛 ) wasconstructed by arranging the vectors in the motion axis
direction of the pin. 𝐴 (3𝑛 3𝑛 ) defines the weight of the
distribution of the acting force for all combinations of pins andmodel nodes. In this study, we assumed that the action of the

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <Figure 3. Example of skin deformation expected by FM1, FM2, DM1, andDM2 algorithms ( 𝑑 2.0 mm, 𝜇 0.6 ). The target deformation wascomputed by the contact simulation.

Figure 4. Example of average deformation and force errors (𝑢 and 𝑓 ,respectively) in progress of time (𝑑 2.0 mm, 𝜇 0.6). The contact state
became stable before approximately 1 s. The difference in the errors
between FM1 and FM2 was relatively small, whereas the difference
between DM1 and DM2 was remarkable.

pins is distributed to the surface of the model by a Gaussian
distribution. On this basis, the Gaussian function values of allneighboring nodes were calculated for each pin, and the weightwas calculated by dividing them by the sum. This method
ensured that the sum of the forces acting on the neighboring
nodes was consistent with that of the forces exerted by the pins.
The distribution of the Gaussian function was empirically set to𝜎 0.8 mm, considering the pin radius (0.5 mm) and internode
distance (approximately 1 mm).
C. Simulation and Matching
The force acting on the finger and deformation of the finger
owing to the contact were calculated using the static linear FEM.
Friction contact with a flat surface was introduced using thepenalty method. Coulomb’s law was assumed for the friction
model. That is, the node that was in contact with the plane waseither in a static friction state or a dynamic friction state. The
transition from the static friction state to the dynamic friction
state occurred when the friction limit was exceeded, and thetransition from the dynamic friction state to the static friction
state occurred because of a stall in the friction speed. The time

6

step of the calculation was 4.17 ms (i.e., 240 steps/s). This timeinterval may not necessarily be sufficient for the accurate
reproduction of the transient state, but it is expected to be
sufficient for the calculation of the friction state following thetransition to a stable dynamic friction.
The device output 𝑓 was calculated using the method
described in Section III for the target states (i.e., 𝑓 and 𝑢) thatwere obtained by the simulation. As an example, thedeformations due to FM1, FM2, DM1, and DM2 at𝑑 2.0 mm
and 𝜇 0.6 are presented in Fig. 3, and their errors are
depicted in Fig. 4. In FM2 and DM2, the optimization problem
with inequality constraints was solved by using the interior
point method. The average errors of the displacement and force
in Fig. 4 were obtained as follows:‖𝑢 𝑢∗ ‖ /𝑛,
20
𝑢
∗
‖
‖𝑓
𝑓 /𝑛.
21
𝑓
The target deformation is constrained by the contact surface,
while the expected deformations by FM1, FM2, DM1, andDM2 can cause deformation that extends beyond the surface
because the constraints by the surface are not imposed on thematching algorithms. The calculation results clarified that thedifference between FM1 and FM2 was relatively small, thedifference between DM1 and DM2 was remarkable, and DM1
exhibited a large error compared to the other methods. In DM1,a solution containing a large negative force was obtained in theunconstrained deformation matching and a large error wascaused by clamping this negative force. In the example
illustrated in Fig. 3, the deformation by DM2 appeared to be theclosest to the target compared to those of the other algorithms.
Based on this result, FM1 and DM2 were used for thecomparison in the experiments. FM2 was not used because thedifference from FM1 was relatively small and it could not be
considered as a generalization of the conventional method.
Moreover, DM1 was considered as unsuitable for comparison
because of its large error.
D. Calibration of Elasticity ConstantIndividual differences in the Young’s modulus of the soft
tissue of the fingers were considered by multiplying the target
force by a constant. In the static deformation calculation by thelinear FEM, the stiffness matrix 𝐾 was a constant matrix
proportional to the Young's modulus. In the contact simulation,
1 10 Pa.
the Young's modulus was assumed to be 𝐸
When the actual Young's modulus of the subject was 𝐸, thestiffness matrix was 𝐸/𝐸 times, and the force generated for thesame deformation was also 𝐸/𝐸 times. To estimate theYoung's modulus of the user, the force against the compression
deformation of the finger was measured and the ratio of thisforce to that of the contact simulation was calculated. This ratio
was expected to approximate 𝐸/𝐸 . The device depicted in Fig.
5 was used for the measurements.
The measurements were performed according to thefollowing procedure. The user maintained their finger lifted sothat the finger touched the push plate lightly. The experimenter
pushed their finger down on the stage and moved it until it
touched the electronic scale. The stage was lowered by 0.5, 1.0,

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <7

1.5, and 2.0 mm from the contact state and the respective loads
were read. Subsequently, α that minimized the sum of squaresof the error between 𝛼 times the load and the force bysimulation was determined. In the experiment, the output force
was calculated using the α value obtained for each subject. Fig.
6 presents the measurement results for each subject. The
obtained value of 𝛼 ranged from 0.61 to 0.90, depending on thesubject.
V. EXPERIMENT
Two experiments were conducted. In Experiment 1, theperception of the tangential and normal forces (forces on the x–
z plane in Fig. 1) when friction was applied to the fingers in
only two directions (left or right) was investigated. In
Experiment 2, the perception of the tangential force (force on
the x–y plane in Fig. 1) owing to friction in any direction of theplane was investigated.
A. Experiment 1: Frictional Force
Subjects were presented the stimuli with a duration of 2 s
calculated by the simulation and matching in section IV.A. The
subjects were instructed to provide feedback on the direction
and intensity of the perceived force in the 𝑥 𝑧 plane. Eachsubject performed 112 trials, which included the following
combinations for the contact simulation conditions: depth 𝑑(0.5, 1.0, 1.5, and 2.0 mm; four ways), coefficient of friction 𝜇(0.0 to 0.6 in 0.1 increments; seven ways), friction direction(right/left; two ways), calculation algorithm (FM1/DM2; twoways), and one iteration. A within-subject design was adopted
in which each subject evaluated all the conditions. The
presentation order was randomized to reduce order effects. Fig.
7 depicts the target deformation 𝑢 and expected deformation 𝑢∗

Figure 7. Examples of target and expected deformations of skin in
Experiment 1. The deformation differed depending on the depth of contact𝑑, friction coefficient 𝜇.

for FM1 and DM2 under certain conditions, all of which werein the state of 2 s from the start of deformation. The subject
provided feedback on the perceived force vector on the 𝑥 𝑧
plane by pointing the mouse on the screen. The criteria for theforce intensity were decided by the subjects. For convenienceof providing feedback, a circle centered on the origin of thevector was displayed on the answering screen.
B. Experiment 2: Direction of Frictional Force
Subjects were presented the stimuli with a duration of 2 s
calculated by the simulation and matching in Section IV.A. The
subjects were instructed to provide feedback on the direction
and intensity of the perceived force in the 𝑥 𝑦 plane. Eachsubject performed 144 trials, which included the following
combination of contact simulation conditions: depth 𝑑 (2.0mm; one way), coefficient of friction 𝜇 (0.5; one way), friction
direction 𝜃 (360° in 10° steps; 36 ways), calculation algorithm(FM1/DM2; two ways), and two iterations. A within-subject
design was adopted in which each subject evaluated all theconditions. The presentation order was randomized to reduceorder effects. Fig. 8 presents the target deformation 𝑢 andexpected deformation 𝑢∗ for FM1 and DM2 under certain
conditions, all of which were in the state of 2 s from the start ofdeformation. The subjects provided feedback on the perceived
force vector on the 𝑥 𝑦 plane by pointing the mouse on the

Figure 5. Homemade device for elasticity measurement. The device wascomposed of a stage and scale. The elasticity constant for each subject was
obtained by measuring the force–displacement relationship and fitting thesimulation data to the relationship.

Figure 6. Force–displacement relationship measured for all subjects. The
slope of the curve differed depending on the subject and the resulting factor
α varied from 0.61 to 0.90.

Figure 8. Examples of target and expected deformations of skin in Experiment
2. The deformation differed depending on the direction of motion 𝜃.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <screen. The force standards were determined by the subjects.
For convenience of providing feedback, a circle centered on theorigin of the vector was displayed on the answering screen.
C. Experimental Procedure
The experimental procedure was approved by the Ethics
Review Board of the University of Electro-Communications(approval number: 20065). The experiment was carried out
according to the following procedure:(1) Experiment explanation:
The experimenter explained the purpose of the experimentand obtained consent for participation in the experiment. The
purpose of the experiment was to investigate the difference inperception owing to the difference in the calculation method(i.e., rendering algorithm) in the presentation of the frictional
force by the pin-array-type force feedback device.
(2) Preliminary questionnaire:
Subjects were asked to provide information regarding their
attributes (i.e., age, gender, and dominant hand), VR experience,and haptic device experience.
(3) Measurement of finger elasticity:
The force correction coefficient 𝛼 was calculated by theprocedure described in Section IV.C and input into theexperimental program.
(4) Explanation of experimental program operation:
The experimenter explained the operation for the experiment
to the subjects. The presentation was started by pressing thespace key. The subjects assessed the direction and intensity of
the force acting on the finger. The subjects responded bydrawing a vector on the screen with the mouse. In case of a
mistake in the response, it was possible to return to the previous
attempt by pressing the “n” key. The force acting on the finger
included both compression and friction components. The
subjects provided answers regarding the intensity of the force
according to their own criteria.
(5) Attachment of device:
The device was attached to the right-hand index finger of thesubject. The device was set so that the positional relationshipbetween the finger and device was as close as possible to that
depicted in Fig. 1.
(6) Practice 1:
Practice was performed for Experiment 1. First, the state ofcontact and friction between the finger and surface was
explained while displaying the target deformation of the finger
on the screen. Thereafter, the presentation and response wererepeated until the subject became accustomed to thepresentation and established the criteria for their response, or
until the subject was satisfied. No feedback on the correct
answer regarding the force direction was provided at this stage.
No time limit was imposed for the practice.
(7) Experiment 1:
The subjects wore earplugs and noise-canceling headphones.
All controls for the progress of the experiment were performedby the subjects, except when the subject required help. A break
was provided during the experiment.
(8) Practice 2:
Practice for Experiment 2 was performed in the same manner8

as described in (6).
(9) Experiment 2:Experiment 2 was performed in the same manner as
described in (7).
(10) Post-questionnaire:
The subjects were asked to describe their impressions of theexperiment freely.
VI. RESULTS AND DISCUSSION
The results of Experiments 1 and 2 are discussed in thissection. Ten subjects (mean age: 22.7, all male, all righthanded) participated in the experiment.
A. Experiment 1: Frictional ForceThe results of Experiment 1 are discussed separately for theperception of the friction coefficient and the perception of theintensity. For the perception of the friction coefficient, first, thetangent/normal ratio 𝑘 was obtained from the force vector
answered by the subject. If correct perception occurred, 𝑘 was
expected to match the coefficient of friction 𝜇. Fig. 9 presents
the relationship between the coefficient of friction 𝜇 andperceived 𝑘. In this figure, the perceived 𝑘 is plotted for eachcondition of depth 𝑑 and friction direction (right/left). The
tangent/normal ratio calculated from the x–z component of thesum of the forces of the target is also indicated by the dashedblack line, which is labeled as “target.” Furthermore, thetangent/normal ratios calculated from the x–z components ofthe sums of the forces presented in FM1 and DM2 are indicatedby the dashed blue and red lines, respectively, which are labeled
as “actual.” In DM2, the ratio 𝑘 tended to increase
proportionally to the presented coefficient, whereas in FM1, it
did not. The slope 𝑏 of the regression line was tested for eachcondition (see 𝑡 and 𝑝 in Table I). In DM2, the difference
between 𝑏 and 0 was significant, whereas in FM1, no
significant difference was observed. To confirm the difference
between the algorithms, Welch's test was performed on theslopes of FM1 and DM2 under each condition (see 𝑡 and 𝑝 inTable I). The difference in the slope was significant under all
conditions. The expected slope 𝑏 was 1 when the friction
coefficient was correctly perceived; however, DM2 exhibited amuch larger slope. This is discussed later in this paper.
Fig. 10 presents the results of the perception of the force
intensity. The intensity calculated from the x–z component of
the sum of the target forces is indicated with a black dashed line,
which is labeled as “target.” Furthermore, the intensities
calculated from the x–z component of each sum of the forces
presented in FM1 and DM2 are indicated by dashed blue andred lines, respectively, which are labeled as “actual.” Because
the subjects decided the scale of the vector when answering, it
was not possible to evaluate the perception of the absolute force
intensity. Therefore, the vector length was normalized by
dividing it by the average of the results obtained under all theconditions by a coefficient of friction 𝜇 0 for each subject.
The mean and standard deviations of the plot were calculated
for this normalized intensity. Therefore, this plot was obtainedfor observing the qualitative tendency of the presented force
and perceived force with the change in the coefficient of friction

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <Fig. 9. Relationship between friction coefficient 𝜇 and tangent/normal ratio𝑘. In DM2, 𝑘 tended to increase proportionally to 𝜇, whereas in FM1, it didnot.
TABLE ISTATISTICAL ANALYSIS OF SLOPE OF 𝑘/𝜇. SIGNIFICANT DIFFERENCES IN
THE SLOPES OF FM1 AND DM2 WERE OBSERVED.
right𝒅 = 0.5 mm
𝒃

𝝈

left

FM1

DM2

FM1

DM2

0.209 0.644

3.496 0.917-0.698 0.869

3.88 1.206

𝒕

0.324

3.813

-0.804

3.217

𝒑

0.747

<0.001**

0.424

0.002**

𝒗

122.0

123.6

𝒕𝟏

-14.6

-15.3

𝒑𝟏

<0.001**

𝒅 = 1.0 mm𝒃

𝝈<0.001**

FM1

DM2

FM1

DM2

-0.532 1.117

4.734 4.533-0.752 1.168

5.292 5.691

𝒕-0.842

3.721-1.164

3.712

𝒑

0.403

<0.001**

0.248<0.001**

𝒗

99.6

94.8

𝒕𝟏

-18.4

-19.2

𝒑𝟏

<0.001**

𝒅 = 1.5 mm𝒃

𝝈<0.001**

FM1

DM2

FM1

DM2

0.436 0.575

1.821 2.153-0.174 1.298

3.735 6.215

𝒕

0.963

2.077

-0.255

2.507

𝒑

0.339

0.042*

0.799

0.015*

𝒗

101.9

95.2

𝒕𝟏

-7.0

-11.8

𝒑𝟏

<0.001**

𝒅 = 2.0 mm𝒃

𝝈<0.001**

FM1

DM2

FM1

DM2

-0.126 1.489

3.084 4.362-0.352 1.11

3.922 10.549

𝒕-0.1732.471-0.559

2.021

𝒑

0.863

0.016*

0.578

0.047*

𝒗

109.6

82.1

𝒕𝟏

-11

-10.4

𝒑𝟏

<0.001**<0.001**

9

Fig. 10. Relationship between friction coefficient 𝜇 and normalized
intensity of perceived force. In DM2, the actual force was considerablylarger than the target, but the perceived force appeared not to be particularly
large.

𝜇 . In both the FM1 and DM2 algorithms, there was noconstraint to match the sum of the forces with the target, andthe value could differ from the target. In DM2, the acting force
was considerably larger than the target in terms of the sum of
the acting forces; however, the force was not perceived to beparticularly large by the subject.
B. Experiment 2: Direction of Frictional ForceThe results of Experiment 2 are discussed separately for theperception of the friction direction and the perception of theintensity. First, for the friction direction, Fig. 11 presents a plotof the perceived direction with respect to the presented direction
𝜃. In this figure, the direction of the force calculated in contact
with the target in the x–y plane is indicated by a black dashedline, which is labeled as “target.” The directions calculated fromthe sum of the forces determined by FM1 and DM2 are alsoplotted with dashed blue and red lines, respectively, which arelabeled as “actual.” In both FM1 and DM2, a region existed
around 270° in which the direction could not be correctly
expressed as the total force vector owing to the orientation of
the device pin and the posture of wearing it on the finger.
Because FM1 and DM2 find the optimum output under thisconstraint, the resulting force does not have the component inthis direction, and the direction of total force (i.e. actual) in
Figure 12 is discontinuous around 270°. For each condition of
FM1 and DM2, there was no error bias across all presentation
angles (FM1: 𝑡 0.315, 𝑝 0.753; DM2: 𝑡 -0.006, 𝑝
0.995). However, several significant deviations from the target
angle depending on the azimuth were observed, which areindicated by asterisks in Fig. 11. When comparing the variance> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <10

the absolute force intensity. Therefore, the force intensity was
normalized by division by the average of the results for all angle
conditions for each subject. The mean and standard deviationof the plot were calculated for this normalized intensity. In
DM2, the change in the presentation force with the angle wasrelatively large, which appeared to affect the perception of thesubjects. A significant correlation between the two was
observed for DM2 (𝑟 0.585, 𝑡 4.213, 𝑝 0.01).

Fig. 11. Relationship between presented and perceived angles of frictionalmotion. Significant deviations from the target are marked with asterisks.

Fig. 12. RMS error of perceived direction and deformation/force errors
depending on target angle. No significant correlation was observed, except
for a weak correlation between the RMS error and 𝑓 in FM1.

Fig. 13. Normalized intensity of perceived force depending on target angle.
A significant correlation between the actual and perceived forces wasobserved in DM2.

of the error across all angles between FM1 and DM2, it wasobserved that the variance of FM1 was significantly larger thanthat of DM2 (𝐹 1.348, 𝑝 0.0024). To clarify the effect ofthe matching error on the perception, the root-mean-square(RMS) error in the friction direction and the matching error due
to FM1 and DM2 ( 𝑢 , 𝑓 ) are depicted in Fig. 12. A
correlation analysis was performed for each RMS error andmatching error. No significant correlation was observed, except
for a weak correlation with the force matching error 𝑓 forFM1 (𝑟
0.371, 𝑡 2.330, 𝑝 0.026).
Fig. 13 depicts the results on the perception of the force
intensity. The plot indicates the intensity calculated from the x–
y component of the sum of the target forces with a black dashedline, which is labeled as “target.” Furthermore, the intensity
calculated from the x–y components of the sums of the forces
presented in FM1 and DM2 are indicated by the dashed blue
and red lines, respectively, which are labeled as “actual.” As in
Experiment 1, because the criteria for answering were decided
by the subjects, it was not possible to evaluate the perception of

C. Discussion
In Experiment 1, it was found that the friction coefficient wasdifficult to perceive with FM1, whereas this was possible with
DM2. However, the recognized tangent/normal ratio 𝑘 wasmuch larger than the presented coefficient of friction 𝜇 .
Although the reason cannot be clarified from this experiment,
possible factors include problems with the algorithm, device
characteristics, human perceptual characteristics, and theinfluence of the response method in the experiment. Further
investigation is required in future experiments. Anotherconcern regarding this result is that in DM2, the perceptual
disturbance appeared to increase near the friction coefficient
𝜇0.6. As the friction coefficient increased, the skinunderwent greater shear deformation. One possible reason forthe disturbance is the difficulty of causing large deformations
by the device used in the experiment. As stated previously, theorientation of the device pins was almost perpendicular to theskin surface, and a small tangential force was applied.
Intuitively, it is difficult and inefficient for shear deformation
to be caused by the distribution of normal forces. This notion
may lead us to designing a device with different pin
configurations that is optimal for the DM algorithm. Regarding
the perception of the force intensity, the output force of DM2
was considerably different from that of the target, but the
subject perception was not significantly affected. The results
suggest that tactile sensation, without the sensation of force thatis derived from bathyesthesia, does not provide accurate
information of the force as a sum on the skin surface.
In Experiment 2, it was revealed that the perception of thefriction direction was possible with both FM1 and DM2 to acertain extent, but DM2 provided better perception accuracy
than FM1. Both FM1 and DM2 exhibited large variability nearazimuths of 90 and 270°. This is considered to be due to the
influence of the device pin orientation. As described in Section
IV.B, the pins were arranged on the cylindrical surface and the
direction of the pin operation was the normal direction of thecylindrical surface. Because the device was mounted at an angle
to the friction surface, the component in the 90° (y+) direction
could be expressed even though it was small. However, thecomponent in the direction around 270° (y-) could not be
expressed in principle. This is considered to have made theperception at this angle particularly unstable.
Interestingly, this direction was also recognized in FM1. The
result that the direction could be perceived by FM1 is intuitively
inconsistent with the result that it was difficult to recognize thetangent/normal ratio in Experiment 1. Further investigation isrequired to elucidate the cause, but one possibility is that theclues to which the subject paid attention changed owing to the

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <difference in the tasks. As the task of Experiment 2 included theinformational element of direction selection, it is possible thatthe subject considered mapping from the stimulus to thedirection by using a method that is not simply physical.
VII. CONCLUSION
We have presented a rendering method that considers theconstraints for multipoint tactile displays. We proposed andformulated a deformation matching method that calculates theoutput force to minimize the error between the target skindeformation and skin deformation owing to the device action.
As a generalization of the conventional method, we also
formulated a force matching method that minimizes the force
error.
To compare the perceptual characteristics of subjects by
deformation matching and force matching, surveys of theperception of the friction coefficient and of the perception of
the friction direction were performed in Experiments 1 and 2,
respectively. The results demonstrated that the perception of thefriction coefficient was more sensitive in DM2 than in FM1 andthat the recognition of the friction direction was more accurate
in DM2 than in FM1. The effectiveness of deformation
matching in the friction expression was confirmed based onthese results.
The experimental results are considered to be meaningful interms of the comparison between FM1 and DM2, but further
studies are required to elucidate the human perceptual
characteristics. The reason that the perception of thetangent/normal ratio was much larger than that of the friction
coefficient in the simulation needs to be clarified through anevaluation of the perceptual characteristics using real objects
and experiments, with devices that have different DoF.
A linear FEM was used to model the finger as the first stepof the study. In reality, the skin tissue is nonlinear andnonuniform. Moreover, in the model of the action by the pin, it
was assumed that the pin and skin surface did not slip, but inreality, slipping may occur and the point of action on the skincan be changed. Furthermore, the interference between the baseand the finger and contact between the fastener and the finger
were not considered. It is expected that a more accurate tactile
presentation will be possible by using a model that considers
these factors.
In this study, we used an approach to minimize thedisplacement error. From the perspective of considering thehuman stimulus perception mechanism, as opposed to
optimization by deformation, optimization by strain energy
density weighted on the region where the receptors are
distributed will be an interesting topic. In such a formulation, it
will be necessary to discuss whether the convergence of thesolution search is guaranteed mathematically.
Real-time processing of the calculations was not included inthis study. It is quite difficult to solve nonlinear programming
problems in real time, but we are beginning to study theacceleration of processing using a GPU. Moreover, we are
optimistic about future advances in computers and consider that
optimization problems, for example, will be calculated
efficiently by quantum computers.

11

Within the scope of this study, the experimental task was
limited to the perception of friction. A further advantage of theuse of multi-DoF haptic devices is expected to be therepresentation of complex shapes. Thus, it is necessary toevaluate the effect of the proposed method on shape recognition
tasks.
REFERENCES
[1]

[2]

[3]
[4]
[5]

[6]
[7]
[8]

[9]
[10]
[11]

[12]
[13]
[14]

[15]

[16]
[17]
[18]
[19]

I. Sarakoglou, N. Garcia-Hernandez, N. G. Tsagarakis, and D. G.
Caldwell, “A high performance tactile feedback display and its
integration in teleoperation,” IEEE Trans. Haptics, vol. 5, pp. 252–263,
May 2012, DOI: 10.1109/toh.2012.20.
Y. Ujitoko, T. Taniguchi, S. Sakurai, and K. Hirota, “Development offinger-mounted high-density pin-array haptic display,” IEEE Access, vol.
8,
pp.
145107–145114,
Aug.
2020,
DOI:
https://doi.org/10.1109/ACCESS.2020.3015058.
S.-C. Kim, C.-H. Kim, G.-H. Yang, T.-H. Yang, B.-K. Han, S.-C. Kang,and D.-S. Kwon, “Small and lightweight tactile display(SaLT) and its
application,” in Proc. WHC, Salt Lake City, UT, USA, 2009, pp. 69–74.
K.-U. Kyung and J.-Y. Lee, “Ubi-pen: A haptic interface with texture
and vibrotactile display,” IEEE Comput. Graph. Appl., vol. 29, no. 1, pp.
56–64, Jan.–Feb. 2009.
A. G. Perez, D. Lobo, F. Chinello, G. Cirio, M. Malvezzi, J. S. Martín,
D. Prattichizzo, and M. A. Otaduy, “Optimization-based wearable tactile
rendering,” IEEE Trans. Haptics, vol. 10, no. 2, pp. 254–264, Apr.–Jun.
2017.
M. Verschoor, D. Casas, and M. A. Otaduy, “Tactile rendering based onskin stress optimization,” ACM Trans. Graph., vol. 39, no. 4, Jul. 2020,
Art. no. 90.
A. G. Perez, G. Cirio, F. Hernandez, C. Garre, and M. A. Otaduy, “Strain
limiting for soft finger contact simulation,” in Proc. WHC, Daejeon,
Korea (South), 2013, pp. 79–84.
T. H. Massie and J. K. Salisbury, “The PHANTOM haptic interface: A
device for probing virtual objects,” in Proc. ASME Winter Annual
Meeting, Symposium on Haptic Interfaces for Virtual Environment andTeleoperator Systems, Chicago, IL, USA, 1994, pp. 295–300.
K. Minamizawa, H. Kajimoto, N. Kawakami, and S. Tachi, “A wearable
haptic display to present the gravity sensation - preliminary observations
and device design,” in Proc. WHC, Tsukuba, Japan, 2007, pp. 133–138.
M. Konyo, S. Tadokoro, and T. Takamori, “Artificial tactile feel display
using soft gel actuators,” in Proc. IEEE ICRA, San Francisco, CA, USA,
2000, pp. 3416–3421, DOI: 10.1109/ROBOT.2000.845250.
T.-H. Yang, S.-Y. Kim, C. H. Kim, D.-S. Kwon, and W. J. Book,“Development of a miniature pin-array tactile module using elastic andelectromagnetic force for mobile devices,” in Proc. WHC, Salt Lake City,
UT, USA, 2009, pp. 13–17, DOI: 10.1109/WHC.2009.4810818.
T. Nojima, Y. Ooide, and H. Kawaguchi, “Hairlytop interface: Aninteractive surface display comprised of hair-like soft actuators,” in Proc.
WHC, Daejeon, Korea (South), 2013, pp. 431–435.
A. Frisoli, M. Solazzi, F. Salsedo, and M. Bergamasco, “A fingertip
haptic display for improving curvature discrimination,” Presence, vol.
17, no. 6, pp. 550–561, Dec. 2008.
A. L. Guinan, M. N. Montandon, A. J. Doxon, and W. R. Provancher,“Discrimination thresholds for communicating rotational inertia andtorque using differential skin stretch feedback in virtual environments,”
in Proc. IEEE Haptics Symp., Houston, TX, USA, 2014, pp. 277–282.
Z. F. Quek, S. B. Schorr, I. Nisky, W. R. Provancher, and A. M. Okamura,“Sensory substitution using 3-degree-of-freedom tangential and normalskin deformation feedback,” in Proc. IEEE Haptics Symp., Houston, TX,
USA, 2014, pp. 27–33.
S. Yazdian, A. J. Doxon, D. E. Johnson, H. Z. Tan, and W. R. Provancher,“Compliance display using a tilting-plate tactile feedback device,” inProc. IEEE Haptics Symp., Houston, TX, USA, 2014, pp. 13–18.
J. C. Bliss, M. H. Katcher, C. H. Rogers, and R. P. Shepard, “Optical-totactile image conversion for the blind,” IEEE Trans. Man Mach. Syst.,
vol. 11, no. 1, pp. 58–65, Mar. 1970.
Y. Ikei, K. Wakamatsu, and S. Fukuda, “Vibratory tactile display ofimage-based textures,” IEEE Comput. Graph. Appl., vol. 17, no. 6, pp.
53–61, Nov./Dec. 1997.
H. Kajimoto, N. Kawakami, T. Maeda, and S. Tachi, “Tactile feeling
display using functional electrical stimulation,” in Proc. ICAT, 1999, pp.
107–114.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <[20] K. Salisbury, D. Brock, T. Massie, N. Swarup, and C. Zilles, “Haptic
rendering: Programming touch interaction with virtual objects,” in Proc.
Symposium on Interactive 3D Graphics, Monterey, CA, USA, 1995, pp.
123–130.
[21] C. B. Zilles and J. K. Salisbury, “A constraint-based god object method
for haptic display,” in Proc. IEEE IROS, Pittsburgh, PA, USA, 1995, pp.
145–151.
[22] H. B. Morgenbesser and M. A. Srinivasan, “Force shading for haptic
shape perception,” in Proc. ASME Dynamic Systems and Control
Division, Atlanta, GA, USA, 1996, vol. 58, pp. 407–412.
[23] M. Minsky, O. Ming, O. Steele, F. P. Brooks, and M. Behensky, “Feelingand seeing: issues in force display,” in Proc. Symposium on Interactive
3D Graphics, Snowbird, UT, USA, 1990, pp. 235–241.
[24] F. Barbagli, A. Frisoli, K. Salisbury, and M. Bergamasco, “Simulating
human fingers: a soft finger proxy model and algorithm,” in Proc.
HAPTICS, Chicago, IL, USA, 2004, pp. 9–17.
[25] M. Ciocarlie, C. Lackner, and P. Allen, “Soft finger model with adaptive
contact geometry for grasping and manipulation tasks,” in Proc. WHC,
Tsukuba, Japan, 2007, pp. 219–224.
[26] C. Duriez, H. Courtecuisse, J.-P. de la Plata Alcalde, and P. -J.
Bensoussan, “Contact skinning,” in Proc. Eurographics, Crete, Greece,
2008, pp. 313–320.
[27] A. Talvas, M. Marchal, C. Duriez, and M. A. Otaduy, “Aggregate
constraints for virtual manipulation with soft fingers,” Trans. Vis.
Comput. Graph., vol. 21, no. 4, pp. 452–461, Apr. 2015.
[28] J. Jacobs and B. Froehlich, “A soft hand model for physically-based
manipulation of virtual objects,” in Proc. IEEE VR, Singapore, 2011, pp.
11–18.
[29] C. Garre, F. Hernandez, A. Gracia, and M. A. Otaduy, “Interactive
simulation of a deformable hand for haptic rendering,” in Proc. WHC,
Istanbul, Turkey, 2011, pp. 239–244.
[30] K. Hirota and K. Tagawa, “Interaction with virtual object using
deformable hand,” in Proc. IEEE VR, Greenville, SC, USA, 2016, pp.
49–56.
[31] K. Sato, H. Nii, N. Kawakami, and S. Tachi, “Electrotactile stimulation
based on strain energy density of the fingertip,” in Proc. EuroHaptics,
Madrid, Spain, 2008, pp. 724–729.
[32] D. Maemori, L. B. Porquis, M. Konyo, and S. Tadokoro, "A Multi-DOF
Haptic Representation Using Suction Pressure Stimuli on Finger Pads,"
in Proc. Eurographics, Crete, Greece, 2014, pp. 285-294.
[33] H. Barreiro, S. Sinclair, and M. A. Otaduy, "Ultrasound Rendering of
Tactile Interaction with Fluids," In Proc. World Haptics Conf., Tokyo,
JAPAN, 2019, pp. 521-526.
[34] H. Barreiro, S. Sinclair, and M. A. Otaduy, "Path Routing Optimizationfor STM Ultrasound Rendering," IEEE Trans. Haptics, vol. 13, no. 1, pp.
45-51, Jan.-March 2020.
[35] S. Weinstein, “Intensive and extensive aspects of tactile sensitivity as afunction of body part, sex, and laterality,” in Proc. First International
Symposium on the Skin Senses, Tallahassee, FL, USA, 1968.
[36] J. M. Loomis and C. C. Collins, “Sensitivity to shifts of a point stimulus:
An instance of tactile hyperacuity,” Percept. Psychophys., vol. 24, no. 6,
pp. 487–492, Dec. 1978.
[37] M. Nakatani, H. Kajimoto, N. Kawakami, and S. Tachi, “Tactile
sensation with high-density pin-matrix,” in Proc. APGV, A Coruña,
Spain, 2005, pp. 169–169.
[38] J. H. Kirman, “Tactile apparent movement: The effects of interstimulus
onset interval and stimulus duration,” Percept. Psychophys., vol. 15, no.
1–6, Jan. 1974, DOI: https://doi.org/10.3758/BF03205819.
[39] M. A. Srinivasan and K. Dandekar, “An investigation of the mechanics
of tactile sense using two-dimensional models of the primate fingertip,”
J. Biomech. Eng., vol. 118, no. 1, pp. 48–55, Feb. 1996.
[40] K. Kobayashi and T. Maeno, “Relationship between the structure offinger tissue and the location of tactile receptors (2nd report, methodof dynamic contact analysis and results for contact between the finger
and plane plate),” Trans. JSME, vol. 64, no. 628, pp. 4798–4805, 1998
(in Japanese).
[41] Y. Makino, N. Asamura, and H. Shinoda, “A whole palm tactile display
using suction pressure,” in Proc. IEEE ICRA, New Orleans, LA, USA,2004, pp. 1524–1529.
[42] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The Finite Element
Method: Its Basis and Fundamentals, 7th ed., Elsevier Ltd., 2013.
[43] N. Yamashita, Nonlinear programming, Application optimization series
No.6, Asakura Publishing Co., Ltd., 2015 (in Japanese).
[44] N. Mitsuhashi, K. Fujieda, T. Tamura, S. Kawamoto, T. Takagi, and K.
Okubo, “BodyParts3D: 3D structure database for anatomical concepts,”

12

Nucleic Acids Res., vol. 37, Database issue, pp. D782–D785, Jan. 2009,
DOI: 10.1093/nar/gkn613.
[45] S. Yoshimura, R. Shioya, H. Noguchi, and T. Miyamura, “Advanced
general-purpose computational mechanics system for large-scale
analysis and design,” J. Computat. Appl. Math., vol. 149, no. 1, pp. 279–
296, Dec. 2002.

Koichi Hirota received his B.S. and Ph.D. degrees from the
University of Tokyo, Tokyo, Japan, in 1988 and 1994,
respectively. He then became an Assistant Professor at
Toyohashi University of Technology, Toyohashi, Japan, in
1995. In 2000, he became an Associate Professor at the
University of Tokyo. He is currently a Professor at the
Department of Informatics, Graduate School of Informatics andEngineering, University of Electro-Communications, Tokyo.
His research interests include haptic devices, haptic rendering,and human–computer interaction.
Yusuke Ujitoko received his B.E. degree in mechanicalengineering and M.A.E. degree in interdisciplinary information
studies from the University of Tokyo, Tokyo, Japan, in 2014
and 2016, respectively. He received his Ph.D. degree from the
University of Electro-Communications, Tokyo, Japan, in 2020.
He was a Researcher at Hitachi, Ltd., Tokyo, by 2020, and iscurrently a Researcher at the NTT Communication Science
Laboratories. His research interests include applied haptic
perception and haptic interfaces.
Sho Sakurai received her B.E. degree in social and information
studies from Gunma University, Maebashi, Japan, in 2007, and
M.A.E. degree in interdisciplinary information studies andPh.D. degree in engineering, from the University of Tokyo,
Tokyo, Japan, in 2010 and 2014, respectively. He is currently a
Project Assistant Professor at the Graduate School ofInformation Systems, the University of ElectroCommunications, Tokyo. Her research interests are multimodal/cross-modal interfaces and human–computer interaction.
She is also active as a manga artist for the introduction of thelatest research on virtual reality and artificial intelligence.
Takuya Nojima received his Ph.D. degree from the Universityof Tokyo, Tokyo, Japan, in 1988. He was a researcher at theJapan Aerospace Exploration Agency in December 2008 and iscurrently an associate professor with the Department ofInformatics, Graduate School of Informatics and Engineering,
University of Electro-Communications. He was a Visiting
Associate Professor at the Massachusetts Institute ofTechnology Media Lab from August 2015 to March 2016. His
research interests include augmented sports, entertainment
computing, and haptic and tangible interfaces.