1. technology and computing

力触覚の提示と計算研究委員会 第 14 回研究会 (2015 年 3 月)
PBD による変形物力覚レンダーリングの改善
Improvement of Deformable objects Haptic Rendering based on Position based dynamics
丁海陽 1) , 三武 裕玄 2) , 長谷川 晶一 3)
Ding Haiyang, Hironori Mitake and Shoichi Hasegawa
1) 東京工業大学 大学院総合理工学研究科 知能システム科学専攻
（〒 226-8503 神奈川県横浜市緑区長津田町 4259 R2-626, teikaiyo@haselab.net）
2) 東京工業大学 精密工学研究所
（〒 226-8503 神奈川県横浜市緑区長津田町 4259 R2-624, mitake@pi.titech.ac.jp）
3) 東京工業大学 精密工学研究所, 独立行政法人 科学技術振興機構
（〒 226-8503 神奈川県横浜市緑区長津田町 4259 R2-624, hase@pi.titech.ac.jp）
概要:
We propose an improvement for stable haptic rendering method for deformable objects. The method
we use for deformation calculation is called Oriented Particle (OP). Since OP is based on position
based dynamics (PBD) simulation, we calculate user-input of virtual coupling by a PBD constraint
in OP. This constraint solves haptic input by displacement to match the PBD feature. We call it
”haptic constraint”. The instability problem occurred in explicit integration haptic simulation is
solved by this method and we achieved a great alleviation of oscillation during the interaction. Our
main contribution is an improved of haptic rendering system even using high spring value or low
model stiffness. This work is achieved by taking advantage of PBD, especially for deformable model
interaction.
キーワード： Haptic rendering, Position based Dynamics, Oriented Particle
1.
INTRODUCTION
Since OP is implemented on Position Based Dynamics
The dynamic simulation of haptic interaction with solid
(PBD), user inputs of the previous step are used to start
or deformable objects has been investigated for more
the calculation of the current step. However, in this gen-
than two decades. Applications in virtual prototyping,
eral PBD approach, haptic user inputs cause oscillations
creative modeling and the medical area proved the value
on the contact position when using low model stiffness or
of haptics in virtual reality systems. In order to build a
a high virtual coupling spring value.
deformable object haptic system of rich reality, our research focuses on the Oriented Particle (OP) method,
We propose a improved stable haptic rendering method
which is faster and more robust than high cost methods
following the PBD concept. This method works in a sim-
such as the general Finite Element Model (FEM).
ilar way to the PBD constraint and can be seen as a PBD
Haptic rendering requires a haptic interface, a virtual
constraint design for haptics. We will describe this in de-
environment and a control law linking them. Generally
tail in section III. Our implementation work is depicted
the device end-point is connected to a virtual environ-
in section IV, while evaluation follows in section V.
ment proxy through the control law, called virtual coupling. The conventional virtual coupling scheme uses a
spring and a damper in parallel between the device end-
2.
RELATED WORK
point and the proxy. Forces and moments resulting from
Over the past decades, much attention has been de-
user inputs are then able to be calculated from the vir-
voted to research on haptic simulation with deformable
tual coupling, and applied to both the virtual model and
objects. Recently, Barbic[1] succeeded in presenting a
the haptic interface.
two object real-time haptic interaction with 1 kHz hap-
∗ Constraints projection for all vertices:
tic rate. He achieved this with a reduced FEM, named
St.Venant-Kirchhoff Deformable Model, by using dimen-
gi = (C1 (pi ), C2 (pi ), ...)
sional model reduction. Galoppo[2] gave a haptic sim-
pi = gi
ulation environment that can handle multiple objects.
– End iteration loop
These objects are simulated by a rotational FEM method
– Update velocity for all vertices:
but it is still inevitable to use the GPU to achieve real-
vi (t1 ) = (gi − xi (t0 ))/∆t
time simulation.
– Apply position for all vertices:
As a geometric deformation calculation, Muller[4] pre-
xi (t1 ) = gi
sented the shape-matching method. Muller[6] extended
his method, introducing OP as a sparse arranged struc-
• End main loop
ture based on shape-matching. Since it is a fast and
robust method, OP is easier to implement and simplifies
the simulation of large models compared to general FEM
without reduction or the support of GPU. Our idea is
taking advantage of these features, as we use OP for our
deformation calculation.
Our research starts from three-degree-of-freedom (3DOF) haptic interaction.
Zilles[8] proposed the god-
object as the first application of Lagrange multipliers for
contact in 3-DOF haptic rendering. The proxy method[7]
published by Ruspini is an extension of god-object by using a spherical shape and an updating algorithm similar
to robotic arm motion. We implement this method in
our system. Other methods, such as the ideal haptic interface point (IHIP)[3], which limits the intersection test
only around neighbor primitives, are not considered at
this time.
3.
tained from the external force fi . Then the predicted
position pi is corrected in the constraint projection step
such that it satisfies a set of constraints (C1 (pi ), C2 (pi ), ...),
such as collision or deformation. This step is looped one
or more times. Finally, the corrected predicted position,
called goal position gi , is applied to the vertex position xi
to be displayed while the velocity vi is updated from the
position difference. By this way of simulation, vertex i is
directly projected to a valid location, under the control
of constraints.
3.2
Haptic constraint
As we described before, haptic rendering needs virtual
coupling to link a dynamic model of the virtual proxy
with a haptic interface. Conventionally, the relative mo-
METHOD
In this section, we will begin with a short explanation
of PBD. Then we will introduce our new method.
3.1
First, the vertex position xi is predicted according to
the velocity vi , computed from the the acceleration ob-
Basic position based dynamics
Since any time integration generally computes velocities first from the accelerations and then positions from
velocities, PBD uses accelerations and velocities of the
previous simulation step to predict positions. Then, after revising the positions by constraints, the velocities
of the current step are calculated from the difference of
positions. This is why PBD is called position based.
For the sake of a clear view of PBD flow, the main
PBD simulation loop is as follows:
tion of the virtual proxy and the haptic interface are
linked by a spring and a damper. Then, user inputs are
transformed to forces and moments and applied to the
model. Finally, the output displayed on the haptic interface is equal and opposite to these forces and moments.
In PBD, force based simulation accumulates forces and
moments at the beginning of each simulation loop. That
means user inputs applied to model are generally involved
in the position prediction step. Since they are calculated
from the previous simulation loop, the effect of user inputs occurs one loop later than the other constraints.
This explicit way of integration becomes unstable when
the user input introduces a relatively large effect compared to the effect of PBD constraints. Examples that
can lead to such behavior are the use of a high virtual
• Initialize xi ,vi for all vertices.
coupling spring value or low stiffness deformable models
• Main Loop:
without damping.
– Compute velocity for all vertices:
vi (t0 ) = vi (t−1 ) + ∆tfi
– Predict position for all vertices:
pi = xi (t0 ) + ∆tvi (t0 )
– Iteration Loop:
By following the concept of PBD, we design a new virtual coupling as a constraint to build a non-explicit integration. This constraint is also a part of the constraint
projection loop. Here, user inputs are introduced to the
PBD solver as a spring constraint and projected together
with other constraints. This concept was first proposed
equation corrects the vertex position to the goal position
of the haptic constraint.
Since the PBD constraint is iterated multiple times,
the effect of kc is non-linear. Here we follow the stiffness
transformation equation[6] of PBD to linearize the effect
of kc by kc′ = 1 − (1 − kc )1/ns . Here ns stands for the
iteration time. Therefore we are able to use kc in a linear
form similar to a conventional spring.
As in OP, vertices are controlled by particles and those
particles are corrected in the PBD constraint loop. We
realize the push by applying a weighted distance value
to surrounding particles resulting in a peak shaped force
Figure1: Calculation of haptic constraint: Dis-
effect. Since replacing a force by a distance applied on
tance from proxy position to haptic device end-
particles will cause both translation and rotation in dy-
point is projected to particle translation and ro-
namics, we compute both for each particle, as shown in
tation to compute the goal position of haptic con-
figure 1. They can be formulated as follows:
straint
The push direction n is computed from the proxy position l and the haptic device end-point u:
by Anitescu[9] in a spring model for implicit integration
n=
of dynamic simulation. In PBD, a stiffness coefficient
is used to replace the effect of the spring. We will give
u−l
,
d
(2)
Here, d is the total pushed distance,
an explanation on how to calculate this constraint in the
d = ||u − l||.
next section.
On the other hand, the feedback force to the user is
calculated in the same way as the conventional method.
We start the deduce from the force integration scheme of
dynamics equations:
∫
The same spring coefficient as the spring constraint is
mv =
used here. Therefore the haptic simulation scheme is no
f dx,
Iω = r ×
by explicit integration are eliminated. This is the main
tic constraint”.
4.
IMPLEMENTATION
In this section, we will describe how our method is
(4)
∫
longer explicit, and the unstable interactions generated
advantage of our proposal. We named this method ”hap-
(3)
f dx.
(5)
Besides equation 3 the pushing distance of proxy d can
be also defined for a single particle as:
∫
d =
V · ndx.
∫
=
(v + ω × r) · ndx.
(6)
(7)
implemented in a PBD system especially for OP configV is the total velocity of this movement. Then, by sub-
urations.
stituting the equation 4 and equation 5, the following
4.1
deduction is:
Calculation of Haptic Constraint
The haptic constraint is used to apply haptic external
∫
effects on a virtual model. Specifically for PBD, this
d
=
constraint works by pushing the model vertices forward.
In a simple PBD vertex case, if the proxy tries to push
=
the vertex, the pushed distance is from proxy position l
=
to the location of the haptic device end-point u together
with a stiffness kc like this:
gi = pi + kc (u − l),
So,
and pi is the predicted position of the same vertex. This
∫∫
f dx =
(1)
where gi is the goal position of the vertex in general PBD
∫
∫
1
f dx + (I −1 r × f dx) × r) · ndx (8)
m
∫
∫
∫
n
(
f − I −1 r × n × r · f ) · ndx
(9)
m
∫∫
n
( − I −1 r × n × r) · n ·
f dx
(10)
m
(
d
1
(m
n + I −1 r × n × r) · n
(11)
By integrating the equation 4 and the equation 5 by time
∫∫
we get
mt = n ·
f dx
(12)
∫∫
Ia = r × n ·
f dx
(13)
∫∫
Use equation 11 to substitute
f dx in the equation
above, we get the final equation for a single particle.
t=
1
d
n 1
,
m (m
n + I −1 r × n × r) · n
a = I −1 r × n
d
.
1
(m
n + I −1 r × n × r) · n
(14)
(15)
Since the pushed distance is weighted for each particle
according to the length from proxy to the particle center,
we use these equation in our implemented system:
tp =
1
wp d
n 1
,
m (m
n + I −1 r × n × r) · n
(16)
wp d
,
1
(m
n + I −1 r × n × r) · n
(17)
ap = I −1 r × n
Figure2: Evaluation. The force based simulation
with increasing spring value.
Finally, the PBD goal position and orientation of the
haptic constraint are calculated for each particle, by the
following equations:
gp = pp + kc′ tp ,
(18)
Rp = Qp + kc′ Ap ,
(19)
where gp and Rp are the particle goal position and orientation matrix, pp and Qp are the predicted particle
position and orientation matrix, while Ap is the rotation
matrix transformed from the rotation vector ap . Thus
Figure3: Evaluation. The haptic constraint with
we are able to apply our haptic constraint on OP.
After the constraint loop, the final goal position and
increasing spring value.
orientation are calculated and the proxy position is settled down. The settled distance from proxy to device
1. Use a determinate spring value kf for the force
based method;
end-point is then passed to the haptic interface to calcu-
2. Let the proxy push the particle to deform a dis-
late the device force feedback.
tance hf , then record hf ;
5.
STABILITY EVALUATION
This section describes our method evaluation.
3. Adjust the stiffness kc of the haptic constraint and
We
make the same proxy motion to push the particle
will give a stability comparison of our haptic constraint
deform a distance hc .
method to the conventional force based method after the
explanation of how we set-up our experiments.
5.1
Set up
To implement the conventional force based method,
we treat the user input as the external effect in general
PBD flow, which is calculated in the prediction step. For
a valid and direct comparison, we want to identify the
same spring effect for both the force based method spring
and the haptic constraint spring. Since PBD’s spring
coefficient is a constraint gradient, the coefficient is not
the same as the ordinary force based spring. However,
the linearly corrected coefficient kc is proportional to a
force based spring coefficient. Therefore, we made an
identification simulation for both methods as follows:
Thus kc matches the conventional spring value kf while
hf = hc . We do this three times to get the three matching
spring values for stability experiments.
The model we used for evaluation is a simple box model
of size 2cm×2cm×2cm and weight 8g. Particles are randomly arranged on the box, totaling 125 vertices with 25
particles.
5.2
Simulation Result
Our system is able to simulate 1kHz haptic interation
in real-time by implementing a multi-rate method. We
use 0.01 second as a fixed time step similar to the general
OP system. We use Spidar-G6[10] as our haptic device
in both the evaluation section and the next section.
Figure6: Demonstration. Screen-shots of haptic
Figure4: Evaluation. The force based simulation
interaction. Our system is able to handle stable
with decreasing model stiffness.
haptic interaction with large deformation. The
slime-rabbit model is dramatically deformed by
user’s haptic operation.
6.
HAPTIC INTERACTION
We tried our system on a simplified stanford bunny
model. This model consists of 458 vertices and 80 particles. The sub index from 1 to 4 in figure 6 show the
screen-shots of interaction with the deforming bunny model.
The corresponding time frames are illustrated in figure
7 and figure 8. Figure 7 shows the rendering force and
the position of device, proxy and particle. The figure
shows that a stable haptic rendering is achieved. Figure
Figure5: Evaluation. The haptic constraint with
8 shows computation time for each simulation step. To-
decreasing model stiffness.
tal computation times are less than 10ms on our Core
i5 laptop and simulations are performed in real-time. In
The evaluation is done by the procedure: move device
figure 8, ”Haptic” is the computation time for the col-
end-point downwards, touch the model and move back
lision detection and the proxy update. ”Deform” is the
up.
computation time for the simulation of OP. ”Blend” is
In the first simulation, we compare different spring val-
the computation time for the surface position updates
ues for both the haptic constraint method and the force
from OP. ”Total” includes the graphics rendering and all
based method. The three spring values kf are 100 N/m,
other computations.
200 N/m and 300 N/m, and the corresponding kc are
The user feels the bunny’s surface shape and the elas-
0.075, 0.15 and 0.225. The Y coordinates of the feedback
ticity smoothly without any jitter or oscillations. As we
force, the positions of the device end-point, the proxy
do not employ force shading, the user can feel the edges
and the nearby particle are displayed in figure 2 and fig-
of each triangle. While the proxy update generally works
ure 3. Oscillation occurs for high virtual coupling spring
well, extreme deformation exceptionally causes falling of
values in the force based simulation, while our method
the proxy through the model surface.
stays stable.
In the second simulation, we compare the two methods
7.
CONCLUSION
with decreasing model stiffness values from 0.5 to 0.3
In this paper, we present an improved stable haptic
with the spring value 300N/m. From figure 4 and figure
rendering method based on OP. This method is a rede-
5 we can see that, when using lower stiffness values, the
veloped proxy method in a PBD simulation system. We
oscillation is larger in the force based method, while our
found that a force based implementation will cause oscil-
method shows steady feedback.
lation during simulation. This problem can be solved by
a new way of force calculation based on PBD constraint.
We realize this by designing a haptic constraint which is
similar to the other constraints in PBD. As a result of
[1] Anitescu M. , Potra A. F.: ”A time-stepping method
for stiff multibody dynamics with contact and friction”, Intl. J. for Numerical Methods in Engineering
2002 55, pp.753-784, 2002.
[2] Barbic J., James D. L.: Six-dof haptic rendering
of contact between geometrically complex reduced
deformable models, IEEE Trans. Haptics, vol. 1, no.
1, pp. 1-14, 2008.
[3] Galoppo N, Tekin S, Otaduy MA et al: Interactive
haptic rendering of high-resolution deformable objects, Virtual Reality (LNCS), 4563:215-223 2007.
[4] HO, C., BASDOGAN, C., AND SRINIVASAN, M.
Figure7:
Demonstration.
The force feedback,
proxy position, device position and nearby particle position during the interaction.
The four
screen-shots are marked.
A.: An efficient haptic rendering technique for displaying 3D polyhedral objects and their surface
details in virtual environments. 1999. Presence in
Teleoper. Virtual Environ. 8, 5, 477 491.
[5] MULLER,
M.,
HEIDELBERGER,
B.,
TESCHNER, M.: Meshless deformations based on
shape matching, In Proc. SIGGRAPH, 471-478
2005.
[6] MULLER, M., HEIDELBERGER, B. J,Hennix.
Ratcliff, Position Based Dynamics, J. Vis. Commun., vol 18, issue 2, pages 109-118, 2007.
[7] MULLER, M., AND CHENTANEZ, N.Y: Solid
Simulation with Oriented Particles, ACM Trans.
Graph,92:1-92:10. 12,13, 2011.
[8] Ruspini D. C., Kolaroy K., Khatib O.: The Haptic
Display of Complex Graphical Environments, COMFigure8: Demonstration. The computation time
for each step during the interaction.
The four
screen-shots are marked.
PUTER GRAPHICS Proceedings 345-352. 1997.
[9] Zilles, C. B., Salisbury, J. K.: A Constraint-based
God-object Method for Haptic Display, ASME Haptic Interfaces for Virtual Environment and Teleop-
our evaluation, the oscillation is greatly alleviated with
erator Systems 1994, Dynamic Systems and Control
large spring values or low model stiffness.
1994 Chicago, Illinois, Nov. 6-11, vol. 1, pp.146-150.
[10] S. Kim, J. Berkley, and M. Sato. A novel seven de-
8.
LIMITATION & FUTURE WORK
Since this work is only focused on three-degree-of-freedom
haptic interaction, further arrangements and variations
are needed for six-degree-of-freedom simulation. Although
our code and system are not yet optimized for practical
application, our research makes stable haptic implementation of the fast and robust OP deformation method
possible. Conventional haptic implemented virtual reality systems are able to use OP as a simple and fast way
to simulate solids. We hope our work can be a reference
for future exploration.
ACKNOWLEDGMENT
This work is support by JSPS KAKENHI Grant Number 26280072.
参考文献
gree of freedom haptic device for engineering design.
VIRTUAL REALITY, 6(4):217 228, 2003.