Interactive Generalized Penetration Depth Computation for Rigid andArticulated Models using Object Norm
Min Tang and Young J. Kim
Dept. of Computer Science and Engineering, Ewha Womans University, Seoul, Republic of Korea
We present a novel, real-time algorithm to accurately approximate the generalized penetration depth (PDg ) between two overlapping rigid or articulated models. Given the high complexity of computing PDg , our algorithm
approximates PDg based on iterative, constrained optimization on the contact space, defined by the overlapping objects. The main ingredient of ouralgorithm is a novel and general formulation of distance metric, the object norm, in a configuration space for articulated models, and a compact
closed-form solution for it. Then, we perform constrained optimization, by
linearizing the contact constraint, and minimizing the object norm undersuch a constraint. In practice, our algorithm can compute locally optimal
PDg for rigid or articulated models consisting of tens of thousands of triangles in tens of milliseconds. We also suggest three applications using PDg
computation: retraction-based motion planning, physically-based animationand data-driven grasping.
Categories and Subject Descriptors: I.2.9 [Artificial Intelligence]: Robotics—kinematics and dynamics; I.3.5 [Computer Graphics]: Computational
Geometry and Object Modeling—physically-based modeling; I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—animation;
I.6.8 [Simulation and Modeling]: Types of Simulation—animation
General Terms: Algorithms, Experimentation, Theory, Verification

The need for PD calculation arises in many graphical applications. For instance, in penalty-based haptics or dynamics simulations, collision-response forces are computed and applied based on
the PD amount between interacting objects [Moore and Wilhelms
1988; Kim et al. 2003]. Furthermore, for constraint-based dynamics, PD can be used for simulation stabilization, to alleviate theconstraint drifting problem [Baumgarte 1972; Hart et al. 2002]. In
retraction-based motion planning, PD is used to solve the very challenging narrow-passage problem, such as the alpha puzzle [Zhang
and Manocha 2008] or to determine the path existence [Zhang et al.
2008].
However, as pointed out by [Zhang et al. 2007a], translational PD
tends to over-estimate the amount of interpenetration, since the separating motion is limited to only translation. This problem becomeseven worse for an articulated model, as illustrated in Fig. 1 (a). In
this example, the red object has a small amount of overlap withthe blue obstacle, but the result of the translational PD (the yellow
object) can be relatively large. However, by introducing a small angular displacement at both ends of the links for the red object, one
can separate it from the obstacle with a small motion, as shown in
Fig. 1 (b).

Additional Key Words and Phrases: Animation, Dynamics, Penetration
Depth, Collision Detection, Articulated models

1.

INTRODUCTION

Measuring the amount of interpenetration between two overlapping
objects is an important problem in computer graphics and interactive techniques. A well-known distance measure to serve such a
purpose is penetration depth (PD), which is often defined as theminimum amount of translation to separate one object from another [Cameron and Culley 1986; Dobkin et al. 1993], also known astranslational PD.
Email: {tangmin,kimy}@ewha.ac.kr∗ Young J. Kim is the corresponding author.
Permission to make digital or hard copies of part or all of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that copies
show this notice on the first page or initial screen of a display along withthe full citation. Copyrights for components of this work owned by others
than ACM must be honored. Abstracting with credit is permitted. To copy
otherwise, to republish, to post on servers, to redistribute to lists, or to use
any component of this work in other works requires prior specific permission and/or a fee. Permissions may be requested from Publications Dept.,
ACM, Inc., 2 Penn Plaza, Suite 701, New York, NY 10121-0701 USA, fax+1 (212) 869-0481, or permissions@acm.org.
c YYYY ACM 0730-0301/YYYY/13-ARTXXX $15.00
DOI:http://dx.doi.org/10.1145/XXXXXXX.YYYYYYY(a) Translational PD(b) Generalized PD

Fig. 1. Different PD Results. (a) A red, movable articulated body with 3
links is collided with a U-shaped, blue obstacle. The yellow object is theresult of the translational PD. (b) The green object is the result of the generalized PD.

As a remedy to this overestimation problem in translational PD,[Zhang et al. 2007a] proposed a new penetration measure, called
generalized PD (in short, PDg ), which is defined as a minimal
rigid motion to separate overlapping objects; here, the minimality is based on some distance metric in configuration space such as
Dg [Zhang et al. 2007a], S and geodesic [Nawratil et al. 2009], displacement [Zhang et al. 2007b; 2007a], and object norm [Zhang
et al. 2007a]. However, the exact computation of PDg may require
an arrangement of high dimensional contact surfaces, resulting inO(n12 ) combinatorial complexity for a rigid model with n triangles
in 3D [Zhang et al. 2007a]. Thus, all existing approaches use anapproximate method, but are still slow for interactive applications.
Finally, no PDg algorithms exist for articulated models.
Main Results In this paper, we present algorithms to accurately approximate penetration depth (PDg ), for rigid and articulatACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

2•

Tang and Kim

ed models consisting of tens of thousands triangles, at interactive
rates. The heart of our algorithm lies in generalizing the definitionof the object norm metric, originally proposed by [Kazerounian and
Rastegar 1992], from rigid to articulated models, and efficiently
calculating it for diverse joint types. Given the high computational
complexity of PDg for rigid and articulated models, our algorithm approximates PDg , by iteratively performing local optimizationon the contact space. We further simplify the optimization step bylinearizing the high dimensional, curved contact space, but employ
novel and efficient contact-space projection based on generalized
continuous collision detection, and seed-point selection based on
successive perturbation and constrained translation, to reduce theapproximation error involved with PDg calculation. As a result, ouralgorithm can compute PDg for rigid and articulated models consisting of hundreds of thousands of triangles in tens of milliseconds. Moreover, we demonstrate the possible application of our PDg
algorithm to motion planning, grasping simulation, and physicallybased animation.

2.

PREVIOUS WORK

We survey existing PD algorithms, and classify them, according to
the underlying measures to define PD.

2.1

Translational Penetration Depth

The translational penetration depth (PDt in short) can be calculated using a Minkowski sum-based formulation, and [Cameron andCulley 1986] is the first algorithm to employ such an approach.
It is known that exact PDt computation requires O(n2 ) and O(n6 )times for convex and non-convex rigid models in 3D, respectively, where n is the number of triangles in the models. [Dobkin
et al. 1993] proposed a more efficient algorithm, using a DobkinKirkpatrick hierarchy-based acceleration structure. A randomized
algorithm was later presented by [Agarwal et al. 2000]. These algorithms are all exact. [Cameron 1997; Bergen 2001; Kim et al.
2004] proposed an approximate algorithm, and [Bergen 2001; Kimet al. 2004] can offer real-time performance. All of these algorithms
can compute PDt only for convex polytopes.
For non-convex models, [Hachenberger 2009] proposed an exact PDt method based on Minkowski sums, but it is relatively slow.
Most practical PDt algorithms are approximate, such as [Kim et al.
2002; Redon and Lin 2006; Lien 2008; 2009] due to the high complexity of computing PDt , but these methods do not run at interactive rates. More recently, [Je et al. 2012] presented a real-timealgorithm based on iterative contact-space projection.

2.2

Generalized Penetration Depth

Very few algorithms exist to compute PDg . [Zhang et al. 2007a]first proposed the notion of PDg , and presented a method to compute only the lower and upper bounds of PDg for rigid models.
Later, the same authors proposed a more efficient technique, using
a contact-space sampling and displacement metric [Zhang et al.
2007a]. [Nawratil et al. 2009] presented methods based on kinematical geometry, using S and geodesic metrics. However, all ofthese methods are rather slow for interactive applications, and it is
not clear whether they are applicable to articulated models. Morerecently, [Pan et al. 2010] used a penetration measure for retractionbased motion planning for articulated models, but this measure is

not a distance metric, and depends on the underlying representation. To the best of our knowledge, no methods exist for PDg computation for articulated models.

2.3

Other Penetration Measures

PDt or PDg are not the only measures to quantify the amount ofinterpenetration for objects. [Ong and Gilbert 1996; Ong 1993] defined growth distance based on the scaling operator, and [Weller
and Zachmann 2009; Allard et al. 2010; Wang et al. 2012] usedpenetration volumes to measure the interpenetration between objects. [Tang et al. 2009b] proposed pointwise PD based on Hausdorff distance for the intersection of surfaces. For non-rigid models,
researchers often rely on distance fields-based representation [Fisher and Lin 2001; Hoff et al. 2002; Sud et al. 2006] for measuringpenetration, and distance fields can be computed rapidly using GPUs. All these measures are generally not directly related to PDt ,nor to PDg , for rigid or articulated models.

3.

GENERALIZED PD COMPUTATION FOR RIGID
MODELS

Intuitively speaking, the generalized PD, PDg , for rigid models canbe defined as the amount of a minimal rigid motion to separate twocolliding objects A and B, under some distance metric Mg , in a
six-degree-of-freedom (6DoF) configuration space. Generally, any
distance metric in configuration space can be used to define PDg .
In our case, we choose object norm σ [Kazerounian and Rastegar
1992] as the underlying distance metric, because the metric has acompact, closed form (see Eq. 1), is invariant with respect to reference frames [Zhang et al. 2007b], and has a physical interpretation for it (see Sec. 3.1), which is useful for such applications as
physically-based animation and motion planning.

3.1

Problem Formulation and Overview

For the rest of the paper, w.l.o.g., let us assume that a movable
polygonal object A is collided with a stationary polygonal object
B in 3D, and the local frame of A is attached to the center of massof A , and aligned with the principal axes of A .
The object norm σA for A placed at q0 and q1 is defined as [Kazerounian and Rastegar 1992; Zhang et al. 2007a]:
σA (q0 , q1 )R
= V1
(x (q0 ) − x (q1 ))2x∈A= V4 Ixx q21 + Iyy q22 + Izz q23 + q24 + q25 + q26(1)

where q = [q1 , q2 , q3 , q4 , q5 , q6 ]T , [q1 , q2 , q3 ] is the vector part of aunit quaternion representing the relative orientation difference between q0 and q1 , [q4 , q5 , q6 ] is the relative position difference between q0 and q1 , x(q) denotes a point on A when placed at q, V isthe volume of A , and Ixx , Iyy and Izz denote the diagonal entries of
the inertia tensor of A .
We can show that σA corresponds to the potential energy E needed
to move A from q0 to q1 , when two points x(q0 ) and x(q1 ) areconnected with a spring with zero equilibrium length and the same
Hooke coefficient k = V2 , for all points x ∈ A :ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

E = 12Rx∈Ak (x (q0 ) − x (q1 ))2 = V1Rx∈A(x (q0 ) − x (q1 ))2 .

(2)

q2 q2

q1

q2

o o

oo o

o

•

Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Normqf qf

Contact Space
Contact Space
Contact Space

o o

LLCS
LLCS
LLCS

q fqf qf
q0 q0 q0

qf

q1

oo

o

o
oo

o

(a) Free-conf. selection

q1 q1
q 0qq01 q 0(b) Contact space proj.

(c) Constrained optimization

q1

q2

3

q2

o

o

o

(d) Re-projectiono(e) Iteration

Fig. 2. Iterative Constrained Optimization for PDg Computation. (a) a collision-free configuration q f is chosen; (b) q f is projected on the contact space toget q0 ; (c) a linearized local contact space is constructed around q0 , and q1 is obtained by constrained optimization on LLCS; (d) q1 is re-projected to the
LLCS LLCS
contact space to yield q2 ; and (e) the algorithm is iterated, until it reaches a locally-optimal qi . PDg ≡ σA (qi , o).

qf

qf

qf

qf

q1 q 0 q1
qq0(5) Iteration: we qiterate0 (3)-(4),0until we cannot further reduce theNote that since we just consider spring potential energy in Eq.2,
the object norm is also the necessary or minimum kinetic energy
object norm σA (Fig. 2 (e)). Then, PDg ≡ σA .
needed to move A from q0 to q1 .
Contact Space
Our iterative algorithm shares some similarity with [Je et al. 2012;
Contact Space
Assuming that the local frame of A initially corresponds to the
Zhang et al. 2007a], but the main contribution of our work is theworld frame o, then PDg using object norm σ can be defined as
constrained optimization that needs to deal with curved contactfollows [Zhang et al. 2007a]:
space unlike [Je et al. 2012] and does not rely on sampling forminimization unlike [Zhang et al. 2007a]. Because of this, our alσ
PDg (A , B)(3)
gorithm is more general and much faster than the existing methods.
= {min {σA (q, o)} kinterior (A (q)) ∩ B = 0,
/ q∈F}oowhere, F denotes the free configuration space. Note that, in general, PDσg (A , B) 6= PDσg (B, A ). Moreover, according to the contact
space realization property [Zhang et al. 2007a], q must lie on thecontact space C of A in the configuration space, and we can showthat PDσg computation can be posed as a minimum distance calculation problem under some non-Euclidean norm k·kM :PDσg (A , B) = min kqkM ,∀q∈C

where, kqkM = qMqT and M is a 6 × 6 SPD matrix,
 4

V Ixx
4


,
V Iyy
M=
4


I
V zz
E

(4)

(5)

where, E is a 3 × 3 identity matrix.
However, since calculating C precisely is prohibitive, requiring
O(n12 ) arrangement computation, we approximate Eq.3 using a linearization of C and performing iterative, constrained optimization
on it, as follows:(1) Free-configuration selection: we select a free configuration q f
close to o (Fig. 2 (a)).
(2) Contact-space projection: we project q f on the contact space
C toward o, and a projected configuration q0 is obtained (Fig.
2 (b)).
(3) Constrained optimization: we build a linearized local contact
space (LLCS) around q0 , and perform constrained optimization on the LLCS to generate a new configuration q1 with alocally-minimal object norm (Fig. 2 (c)). Kinematically, thisstep corresponds to a gliding motion of A on the LLCS fromq0 to q1 while minimizing its spring potential energy, as defined in Eq. 2, when one end of the spring is attached to o.
(4) Re-projection: If q1 ∈/ C due to the linearization and localization, we re-project q1 to C to get q2 (Fig. 2 (d)).

3.2

o

o

o

Algorithm

Now, we explain our PDg algorithm for rigid models in more detail.
Free-configuration Selection In our algorithm, we can use any
method to find a collision-free configuration, potentially close too. As an implementation choice, we use a variant of methods proposed by [Je et al. 2012], such as the centroid difference,
maximally-clear configuration, and motion coherence method. In
particular, motion coherence can be more effectively utilized for
PDg than PDt algorithms by keeping and reusing the orientation
information in the history of motion. Reusing orientation is impossible for PDt algorithms, since they allow only translation as aseparating motion.
Contact-space Projection The objective of this step is finding acontact configuration qc projected from a collision-free configuration q f toward o. Essentially, this step is very similar to what
continuous collision detection (CCD) algorithms do [Zhang et al.
2007b; Tang et al. 2009a]. We use a variant of [Tang et al. 2009a] toperform CCD more efficiently in our case as follows. We first rotate
A (q f ), to have the same orientation like A (o), unless the rotation
creates new collisions, and then perform translational CCD fromq f to o, by ignoring the rotational part of the interpolating motion
in CCD; otherwise, we resort to full CCD. In practice, this approach
is more effective, than executing full CCD from the beginning.
Constrained Optimization The goal of this step is to constructa contact space (LCS) locally around the contact sample, qc , obtained from the projection step, and to find a configuration on LCS that has a minimal object norm. In general, LCS consists of 5dimensional non-linear manifolds that make non-linear constrained
optimization on LCS very challenging for interactive applications.
We circumvent this problem, by linearizing the LCS to LLCS.
We first explain how to build the LLCS, and then show how to perform constrained optimization on it. The non-penetration constraint
in work space around q can be expressed as J0 V ≥ 0, where J0 is theconstraint Jacobian [Catto 2005; Ruspini and Khatib 1997], whose

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

o

4

•

Tang and Kim

row vector of J0 is
J0 i = (ni )T [−r̂i I], i = 1, ..., n(6)where, n is the number of contact points at q, ni is the normal ateach contact point, r̂i is the skew-symmetric matrix representing
a displacement
between the contact point and origin of A , and ωV=is a generalized velocity of A at qc . J0 i is the Plücker
v
coordinate of the line defined by ni and ri .
Then, we use another Jacobian


RJω 0
Jq =
,
0 I

to time. The solution set for Eq.9 is convex and its boundary set
corresponds to the LLCS around the neighborhood of qc . Now, since we define the LLCS in q’s space, the constrained optimization
problem, in terms of object norm, can be posed as projection on aconvex polyhedral cone:(10)

Then, we can find a minimizer, without relying on sampling thecontact space, like [Zhang et al. 2007a]. Since M is positivedefinite, and J has a full row rank after redundant rows are removed,
Eq. 10 can be solved by the Lagrangian multiplier method:dq∗ = M−1 JT λ − qc .

(11)

As shown in [Redon et al. 2002], calculating λ is equivalent to alinear complementarity problem (LCP), formulated as: JM−1 JT λ − Jqc ≥ 0
λ ≥0 
(12)

JM−1 JT λ − Jqc λ = 0.
We use the projected Gauss-Seidel (PGS) method to find a solution
for Eq.12. The PGS method is relatively simple to implement, andfast, compared to other LCP solvers [Jourdan et al. 1998; Catto
2005; Otaduy et al. 2009]. Moreover, since J has a full row rank andM is a symmetric positive definite (SPD) matrix, JM−1 JT is also
an SPD matrix, which makes the Gauss-Seidel method convergent.
Once we have a solution for λ , from Eq. 11, we can get an optimal
configuration on LLCS q∗ = dq∗ + qc .
Note that the orientation component of q∗ is a unit quaternion, andthus (q∗1 )2 +(q∗2 )2 +(q∗3 )2 ≤ 1. If the result of optimization does notsatisfy this, we use a weight 0 < w < 1 to adjust q∗ = wdq∗ + qc , sothat the quaternion component of q∗ becomes a unit vector. Since
our constrained optimization algorithm always keeps the object
norm for dq∗ + qc smaller than that for qc , and also since the object
norm is a convex function, the object norm for q∗ = wdq∗ + qc issmaller than that for qc .

qi

o

o
(a) Collision-free case(b) In-collision caseFig. 3. Re-projection Re-project qi onto the contact space to get a newcontact configuration qi+1 .

4.

(9)

where, J = J0 Jq , and dq is the first derivative of q with respectqi1 qi1

qi1

(7)

To redefine the constraint from work space to configuration space,using V = Jq dq [Mirtich 1996], the constraint becomes

Minimize 21 kqc + dqk2M
Sub ject to : Jdq ≥ 0.

qi

qi1

where R is the rotation matrix of A at qc , and
q20 + q21
q0 q3 + q1 q2 −q0 q2 + q1 q3
2 Jω =−q0 q3 + q1 q2
q20 + q22
q0 q1 + q2 q3  . (8)q0
q0 q2 + q1 q3 −q0 q1 + q2 q3
q20 + q23

Jdq ≥ 0,Re-projection Due to the nature of linearity and localization in
LLCS, the configuration qi resulting from the constrained optimization may not lie on the contact surface,qt 1as illustrated in Fig.
3. There can be two cases for qi : collision-free
and in-collision. In
qt 1contact
qt
the former case, we re-project from qi onto thespace toward o, to get a new contact configuration, qi+1 . In the latter, we
re-project from the last known contact configuration qi−1 towardoqi , to obtain qi+1 .

GENERALIZING OBJECT NORM FOR
ARTICULATED MODEL

In this section, we generalize the definition of object norm for anarticulated model with different joint types, and use it to define andcompute PDg for articulated models in the next section.

4.1

Notations

We introduce some notations used to define object norm for articulated models. We represent an articulated model A made up ofn links with no closed loop. For the sake of simplicity, we assume
that the parent of link Li is Li−1 , even though this can be easily relaxed. The joint that connects two links Li−1 and Li is denoted asJi , and J0 is the joint attached to the world reference frame. For agiven link Li , {i − 1} denotes its local reference frame of the link,and {0} is the world reference frame. For a joint Ji , we call the link
Li and all its descendent links (i.e. Lk , ∀k > i) as the outboard links
Li ≡ {Lk |k ≥ i}, as illustrated in Fig. 4. We refer to the configuration of a link Li as that of its inboard joint Ji .

4.2

Object Norm for Different Joint Types

To begin with, we define the object norm for one joint Ji of thearticulated model, and later extend it to a more general case. As
we change the joint parameter qi of Ji while fixing all other joint
parameters, only its outboard links Li will rigidly move. Thus, wedefine the object norm σi for Ji as the sum of squared distances ofpoints on Li when only qi changes:
σi =

1
V

Z


2
x q0i − x (qi ) dx

(13)

x∈A

1
=
V

Z


2
x q0i − x (qi ) dx

x∈Li

where V is the volume of A , q0i and qi are the parameters for Ji at
two different configurations, and x is a point on Li in workspace.
Now, we will show how to calculate the object norm of Ji when it
is prismatic, revolute, or spherical based on Eq.13.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Norm

Prismatic Joint In this case, since any point in Li travels the same
distance when the joint parameter of Ji changes from qi to q0i , theobject norm for a prismatic joint can be calculated as:
σi =

Vi 0
|q − qi |2 ,
V i

(14)

where Vi is the volume of Li .
Revolute Joint For a revolution joint Ji , we calculate its object
norm in reference frame {i − 1}, since the Euclidean distance isinvariant to the choice of reference frames. Let us say that R is arotation matrix corresponding to the change of the joint parameterof Ji , r is a vector from the origin of {i − 1} to a point in Li , ωiis the rotational axis of Ji in frame {i − 1}, and θ is the rotational
angle of Ji , as illustrated in Fig. 4.

Li 1

Ji

Li

r

i
Fig. 4. Revolute Joint and Notations.

We can calculate the object norm for a revolute joint Ji as follows:
2
Rx(q0i ) − x(qi ) dx
σi = V1
(15)
Rx∈Li
= V1
(Rr − r)T (Rr − r) drr∈Li

By using the Rodrigues’ formula, we get  Rσi = V4 sin2 θ2rT r − rT ωi ωiT r dr
r∈Li

(16)

where, θ = q0i −qi and by setting r = [xk , yk , zk ]T , we finally obtain:  R
σi = V4 sin2 θ2ωiT rT rI − rrT ωi dr
r∈L
i
 2 2
 
yk + zk −xk yk −xk zk
 
R
4
θ
2
T
2
2
 −xk yk x + z −yk zk  dk ωi= V sin 2 ωi kk[xk ,yk ,zk ]T ∈Li−xk zk −yk zk xk2 + y2k
 = V4 sin2

θ2

5

from the base toward the terminal link, and from the parent toward
the child link; i.e. q0 → q1 → · · · → qn−1 . Let us call this type ofparameter change sequential change, and we assume the sequentialchange of parameters is used to define the object norm for articulated models for the rest of paper. Then, the object norm for A undersequential change from q to q0 can be calculated as: n−1
σA q, q0 = ∑ σi .

(19)

i=0

Note that Eq. 19 is a generalized version of Eq. 1, since Eq. 19
can be reduced to Eq. 1 by considering the six-degree-of-freedom
motion of a rigid body as a combination of three prismatic joints
and one spherical joint defined on its local frame.
We can also derive the physical meaning of Eq. 19, as follows. Since, for a single joint, its object norm corresponds to the necessary
potential energy to change its configuration, similarly to the rigid
case, the object norm of the whole articulated model is equivalent to
the necessary potential energy to change its joint parameter undersequential change. The sequential change also penalizes the motionof links that are closer to the base link, than to the terminal. In ouralgorithm, any link can be chosen as a root (base link) that requires
a minimal motion compared to other links. In practice, one may
choose the link as a root that has the largest mass (e.g. the torso
for a human character, and the palm for a hand model), so that theresulting motion would look natural.

i{i  1}•ωiT Ii ωi ,5.

GENERALIZED PD COMPUTATION FOR
ARTICULATED MODELS

To compute PDg for an articulated model, we employ a five-step
approach similar to the one for rigid models, explained in Sec. 3.1.
However, the way that these steps work for articulated models isquite different from rigid models. The main hypothesis underlying
our iterative constrained optimization is the contact-space realization theorem, and we start this section by showing this theorem forarticulated models, and then explain the details of each optimization step.
Theorem 1 For given colliding articulated models A and B
with no closed loops, and A being initially placed at an incollision configuration qin , if q is a configuration of A that realizes
PDg (A , B), then q is on the contact space induced by A and B.
Please refer to Appendix A for proof.

(17)
where, Ii is the inertia tensor of Li in reference frame {i − 1}.
Spherical Joint For a spherical joint, we again use ωi , θ as therotational axis and angle, respectively. However, unlike the revolutecase, ωi is not fixed, and let us denote  the vector part of a unit
quaternion as Q = [q1 , q2 , q3 ]T = sin θ2 ωi . Then, the object norm
for a spherical joint Ji can be written as:
 
θ
4
4
σi = sin2
ωiT Ii ωi = (Q)T Ii Q.
(18)
V
2
V
Full-body Case Now we calculate the object norm of the wholearticulated model A whose joint parameters change from q to q0 ,5.1

Free-configuration Selection

The approximation quality of our iterative, constrained optimization method depends on the choice of a starting, collision-free configuration. One could have used a similar technique like the one
in Sec. 3.2 for a rigid body, but it can be inefficient for articulated
models, so that the chosen, free-configuration may be far from being optimal. Moreover, self-collision between links and joint limit
also need to be taken into account. In this section, we propose threenovel techniques, successive perturbation, constrained translationand motion coherence, to choose a free configuration for an articulated model A that is initially in collision.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

6•Tang and Kim

Successive Perturbation The successive perturbation method tries
to get a free configuration for A , by perturbing a link, as well as its
descendants that are in collision. Moreover, we want to find a free
configuration as close as possible to the given, in-collision configuration, in terms of the object norm value (i.e. smaller object norm
value).
According to our definition of the object norm for an articulated
model, the object norm for a single joint is proportional to the volume or inertia tensor of all its outboard links (i.e. Eqs. 14 and 17).
Then, for the same amount of value changes in a joint parameter,
the object norm for a given joint Ji is always less than that of its ancestor joints (Jk , ∀k < i). Thus, we start perturbation from theterminal (Ln−1 ), supposedly having the smallest object norm value,
and move toward the base link (L0 ), that is in collision.
First, we use a tree T to represent the structure of an articulated
model A , where a node ni in T represents a link Li , and an edge
connecting ni−1 and ni denotes the inboard joint Ji of Li , as illustrated in Fig. 5. In the figure, the red node denotes an in-collisionlink, and the blue one denotes a collision-free link. Our successive
perturbation works as follows:(1) We add all in-collision nodes (e.g. red nodes in Fig. 5 (a)) intoa list L, and sort them in decreasing order of their tree depth.

Algorithm 1 Successive Perturbation
Input: L: a list of collision nodes in decreasing order of depth.
Output: Success or Fail.
1: while (kLk 6= 0) do
2:
n = L. f irst();3:
L.pop f ront();4:
if (!Perturbation(n)) then5:
if (n = root) then6:
return Fail;7:
end if8:n = n.parent();
9:
if (kn.children()k > 2 ) then10:
return Fail;11:
end if12:
if n ∈/ L then
13:
L.push f ront(n);14:
end if15:
end if
16: end while
17: return Success;(2) We start perturbation from the first node in L; e.g. node 5 inFig. 5(a). The perturbation for a node is executed by linearlyvarying its inboard joint parameter within some small range,
and checking for collision for the node (link), as well as for allof its outboard links.
(3) The perturbation for a node is considered successful, when
the link itself, as well as all of its descendent links, become
collision-free; otherwise, we move on to perturbing its parent. In Fig. 5(a), if perturbation for node 5 fails, we will perturb node 4. To prevent the perturbation from affecting othercollision-free descendent nodes, we stop perturbation for a node with more than one child node. An example in Fig. 5 (a) isthat if perturbation on node 4 fails, perturbation on node 3 willcontinue, but if perturbation on 3 still fails, we stop perturbing
node 1, since it may affect the collision-state of node 2, and we
declare that the successive perturbation fails.
(4) The successive perturbation will stop when L becomes empty, or a node with more than one child is encountered. If L
becomes empty, a free configuration is found; otherwise, successive perturbation cannot find a free configuration, and weneed to resort to other methods, such as constrained translation
or motion coherence.
Note that, when the original input of the algorithm contains an incollision node with multiple descendants (e.g. node 1 in Fig. 5(b)), our perturbation algorithm immediately stops (step 4 of theabove procedure). In this case, before executing successive perturbation, we preprocess the input tree structure, by perturbing all thein-collision nodes with multiple descendants in ascending order of
their height; for instance, in Fig. 5 (b), we perturb the in-collisionnodes in the order 0 and 1. If this preprocess succeeds, we moveon to the normal successive perturbation; otherwise, we skip it andresort to other selection methods. Our successive perturbation algorithm is also summarized in Alg. 1.
Constrained Translation Since the successive perturbation
method is a greedy method, it does not always guarantee the generation of a collision-free configuration. In this case, we switch to the

In‐collision
Collision‐free

0
6

1
2

0

7

3

8

4

9

6

1
2

7

3

8

4

9

5

5
(a)

(b)

Fig. 5. Tree Structure of Articulated Model. (a) The in-collision nodes
have a single descendent, or (b) multiple descendants. The latter requires
preprocessing, before executing successive perturbation.

constrained translation method, which is relatively simple to implement. If A is initially self-collision free, then we consider theentire A as a single rigid body, and use the method explained inSec. 3.2. However, if A has a joint constraint in terms of translation for the base link, we use a variant of the centroid-difference
method [Je et al. 2012], by projecting the centroid-difference direction onto the constraint-free direction; for instance, the base ofa puma robot may move around only on the X-Y plane, and X andY will be the constraint-free directions. Thus, the new centroid osof A using constrained translation is calculated as:o − oB(rA + rB ) A· ai aikoA − oB ki∈L0

os = oA + ∑(20)

where, oA and oB are the centroids of A and B, ai is theconstraint-free direction for the base link L0 , and rA and rB are
the diameters of the enclosing spheres for each object. Finally, acollision-free configuration for A is obtained, by translating A
such that its centroid becomes os .

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

•Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Norm

If an articulated body has self-collision initially, and successive perturbation fails, we choose a self-collision free configuration frommotion coherence (e.g. time-dependent physics simulation) or userprovided initial pose (e.g. T-pose in motion capture), and then apply
the constrained translation method.
Motion Coherence If the underlying environment exhibits motion
coherence, we can exploit this to find a collision-free configuration.
This technique was also used for a rigid body in Sec. 3.2, but themain difference here is that we need to cache not only collisionfree, but also self-collision free configurations from the history of
motion.

5.2

Contact-space Projection

L2

L1

L0

L1

L0
obstaclebt lobstacle
bt l

(a)

(b)

Fig. 6. An Example of CCDg . (a) When L1 is in contact, L0 cannot move.
(b) L2 can still move even after L1 is in contact.

linearized. For a spherical joint with a joint limit, it is treated as a
combination of three revolute joints.

The objective of contact-space projection (CSP) is to find a contact configuration qc as close as possible to the given, in-collisionconfiguration qin 1 , in terms of object norm. The CSP for an articulated model is more delicate than for the rigid counterpart, since it
is desired that qc is not only a contact configuration, but should also
have a minimal parameter difference from qin , which will yield aminimal object norm.
Our novel CSP technique consists of two steps: the generalized
continuous collision detection (CCDg ) and the greedy optimization.
These two steps use different types of underlying motion; the former one uses a linearly interpolation motion just like normal CCD
operations [Zhang et al. 2007b], and the latter uses the sequentialchange.
Generalized Continuous Collision Detection We explain how
CCDg works first. We create a motion linearly interpolating between a collision-free configuration q f and in-collision configuration qin in a configuration space. Then, we linearly and simultaneously change the joint parameters, until we encounter new collisions, or the joint parameter reaches the corresponding one in qin .
We refer to this process as generalized continuous collision detection (CCDg , in short), which is a generalized version of continuous
collision detection that calculates only the first time of contact for
the entire object [Zhang et al. 2007b].
More formally, the configuration qi for a link Li obtained by CCDgis:qi = qi,0 + τi dqi , τi ≡ max ∀t ∈ [0, 1] Li qi,0 + tdqi ∈ F
t(21)
where, dqi = qi,1 − qi,0 and qi,0 and qi,1 are the start(collision-free)and goal (in-collision) configurations of Li , respectively. Then, theconfiguration of the entire body after CCDg is q = [q0 q1 · · · qn−1 ].
A simple example is illustrated in Fig. 6. In the figure, the blue
object is placed at start q0 , and the red one is at goal q1 . Note thatthe link L1 is in contact, but if we still keep moving its ancestor link
L0 , a collision can be generated (Fig. 6 (a)); however, we can stillmove its descendent link L2 , until they have the same parameters
as the goal configuration q1 (Fig. 6 (b)).
The CCDg algorithm can be efficiently implemented using a variation of CCD function like [Zhang et al. 2007b]. Here, we give
more details about the interpolating motion used for CCDg . For
revolute and prismatic joints, the joint parameter is linearized, andfor a spherical joint, the orientation component of quaternion is
1 Initially, q = o.
in

L2

7

We maintain a set of links L that contains the links whose joint
parameters can be changed. Initially, L includes all the links in A .
Then, we invoke CCD for a sub-structure of A consisting only of
the links in L, and then find the links that are in contact, and theselinks, as well as their ancestors in the kinematic chain (called N),
are deleted from L, since the parameter change in N may yield
collision. We continue to do this process for the rest of the links in
L, until L becomes empty. The CCDg algorithm is summarized inlines 1-6 in Alg. 2. Here, the CCD function (q, N)=CCD(L, q f , qin )finds a configuration q for the links in L at the first time of contact
when L moves from q f to qin , and also saves the in-contact links
and their ancestors at N. Also note that, during the CCD check, only
the parameters of L change, and the other link parameters remain
the same as the ones in q f .
Greedy Optimization Even after CCDg is finished, it may be still
possible to decrease the object norm of A . For instance, as illustrated in Fig. 7, even after CCDg (the blue object in Fig. 7 (a) is theresult), we can still move the base link, until it contacts the obstacle(the blue object in Fig. 7 (b)), thereby further reducing its object
norm. Thus, greedy optimization is followed after CCDg as a second step of CSP, to reduce the object norm. Let us call the result of
CCDg as q0f , which is a new start configuration for greedy optimization. Then, each link, one by one from base to terminal, will change
its parameter linearly from q0f to qin (i.e. sequential change), until
it creates a collision. The reasoning behind the sequential parameter change for greedy optimization is that, with the same amountof parameter change, a parent link has a higher reduction on object
norm than its children.
In more detail, we use a list T to keep all the links whose parameters
can change, and sort them in increasing order of their depth value;
initially, T includes all the links in A , sorted in increasing orderof depth. Then, we visit each link in T one-by-one, and change its
link parameter while keeping the configurations of other links sameas q f , until we find a configuration q at the first time of contact.
This will be repeated until all links in T are visited. The greedy
optimization is also summarized in lines 7-14 in Alg. 2.
Reverse projection As an articulated model moves from q f to qin ,there can be an obstacle blocking the motion path, which makes
the contact configuration obtained from contact-space projection
far from being optimal, as illustrated in Fig. 8; in the figure, qc
has a smaller object norm than q0c . To deal with this problem, we
may reverse the projection direction from qin to q f , as opposed toq f to qin . Theoretically speaking, performing this reverse projection precisely can be as hard as PDg computation itself; the reverse
projection is also known in the literature [Dobkin et al. 1993] as di-

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

8•Tang and Kim

optimization is performed on it to generate a better contact configuration, in terms of the object norm.
Let us begin by explaining how to set up the local contact-space (orcontact constraint), and linearize it (i.e. LLCS). At a configuration
q with m contact points pi and m contact normals ni in workspace,
we have the following:

Obstacle

Obstacle(a) CCDg(b) Greedy optimization

Fig. 7. A Greedy Algorithm to Get a Contact Configuration. (a) CCDg
detects the first contact (self-collision) between the first link and the lastlink; the contact point is circled in dotted line. (b) By using the greedy
method, the rest of the links can still move until a contact with the obstacle
is found to yield a smaller object norm.

Algorithm 2 Contact-space Projection
Input: L, T, a collision-free configuration q f , and an in-collisionconfiguration qin .
Output: A contact configuration q.
1: {First step: CCDg }2: q = q f ;3: while (kLk 6= 0) do
4:(q, N)=CCD(L,q,qin );5:
L.remove(N);6: end while
7: {Second step: greedy optimization}
8: while (kTk 6= 0) do
9:
l = T. f irst();10:
q0l ≡ ql,in ; {lth parameter of qin }11:
q0i ≡ qi , ∀i 6= l ;12:(q, −)=CCD(l, q, q0 );13:
T.remove(l);14: end while
rectional penetration depth. Thus, we approximate the reverse projection, by discretely sampling the reversed path, and performing
static collision detection along the discretized path. We have foundthat this method generally performs well in practice. Finally, we replace the reverse projection CCDrg with CCDg in lines 4 and 12 in
Alg. 2, to get a better contact configuration.
qf
qc '

qc
qin

C‐Obstacle

Fig. 8. Reverse Projection in Configuration Space. Since qc is more
optimal than q0c , the projection direction is reversed from qin to q f .

5.3

Ji dq = dpi
nTi Ji dq = nTi dpi ≥ 0Constrained Optimization

Similarly to the rigid case, once a contact configuration qc is found,the contact constraint is set up and linearized, and then constrained(22)

In Eq.22, J maps the configuration change to the position change inworkspace. Then, we can get an m× | q | contact matrix C whose
ith row is nTi Ji , and the constraint inequality as follows:
Cdq ≥ 0(23)

Note that Eq.23 is a generalization of Eq. 9 for an articulated model.
Now, our constrained optimization problem can be formulated as:n−1

Minimize σ (q) = V1 ∑ σi ,i=0
Vi |qi |2 , if Ji is a prismatic jointwhere σi = 4 sin2 q2i ωiT Ii ωi , if Ji is a revolute joint
4 (Qi )T Ii Qi , if Ji is a spherical joint
subject to : C(q − qc ) ≥ 0(24)

where qc is the contact configuration obtained from contact-space
projection. To rewrite Eq.32 for a revolute joint in terms of

−1

d 2 sin q2i , we get
2 sin q2i , since dqi = cos q2i
  q −1  q 
i
i
Cdqi = C cos
d 2 sin,(25)
2
2
and create a new variable X for Eq. 32:
qi , ifJi is prismatic
Xi = 2 sin q2i , if Ji is revolute(26) Q , if J is spherical.
iiThen, we rewrite Eq.25 to C0 dX ≥ 0, and Eq.32 to:
Minimize 12 kXk2H
subject to : C0 dX = C0 (X − Xc ) ≥ 0
where H is a symmetric matrix, and Vi , if Ji is prismatic
Hii = ωiT Ii ωi , if Ji is revolute 4I , if J is spherical.
i
i

(27)

(28)

Since H is an SPD matrix, we use the Lagrange multiplier method
to solve Eq.27:
TX∗ = H−1 C0 λ ∗ .
Also, λ can be calculated from an LCP formulation: C0 H−1 C0T λ − C0 Xc ≥ 0
λ ≥0 0 −1 0T
C H C λ − C0 Xc λ = 0.

(29)

(30)

Using the projected Gauss-Seidel method, we can get a solution λ ∗of the above equation.
Kinematic constraints When an articulated model has other kinematic constraints such as joint limit, we need to take this into account in the optimization. For revolute and prismatic joints, the

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Norm

joint limits are represented as angles or displacements that are already included as variables in object norm definition. Thus, these
constraints can be easily handled by truncating the resulting optimized configuration within the joint limit. However, for a sphericaljoint, when the joint limit is represented as an interval using Euler angles, since there is no tight interval-to-interval relationshipbetween quaternion and Euler angles, it is hard to directly incorporate the Euler angle representation into our constrained optimization framework. In this case, we represent the spherical joint as
three revolute joints. Note that the sum of object norms for three
revolute joints is an approximation of that for the spherical joint.

5.4

Re-projection

Similarly to the rigid case, the resulting configuration q of constrained optimization may not lie on the contact space. We need to
re-project q to the contact space. If q is in free space, the contact
space projection is invoked again from q to o; otherwise, from qc toq, where qc is the last known contact configuration obtained fromthe iterative, constrained optimization.

6.

RESULTS AND DISCUSSIONSWe show implementation results of our PDg algorithm, and compare them against other known algorithms. We also demonstratesome applications of our fast PDg algorithms, those of retractionbased motion planning, physically-based animation and datadriven grasp planning. Finally, we discuss the computational complexities, as well as some limitations of our algorithm.

6.1

Implementation and Benchmarks

We implemented and benchmarked our algorithm using C++ on a
Windows 7 PC, equipped with an Intel i5 2.40GHz CPU and 4.0GB memory. We also use the public-domain CD and CCD libraries
CATCH [Zhang et al. 2007b], C2 A [Tang et al. 2009a], and PQP[Larsen et al. 2000] in our implementation, and modify them tosuit our purpose. No parallel processing technique was used, eventhough it may be possible. For free-configuration selection, if motion coherence is present in the simulated environment, it is exploited; otherwise, the centroid-difference method was used for rigid
models, and for articulated models, the successive perturbation is
executed first, and if it fails; the constrained translation is followed.
PDg Computation for Rigid Models The benchmarking models
are shown in Fig. 9, whose complexities range from 20 to 174K
triangles. To measure the performance of our algorithm, we translate A toward B while changing its orientations 200 times, and theaverage performance statistics for some models using two different
free-configuration selection methods (motion coherence and centroid difference) are obtained as shown in Table.I. A typical case ofthis benchmarking scenario is shown in Figs. 10 (a) & (b), where
the red object represents the colliding one, and the green one denotes the result of PDg .
We also compare our method with [Zhang et al. 2007a], which is
known to be the fastest practical PDg algorithm for rigid models, while using the same free-configuration selection scheme. The
comparison results, including the average relative error and average timing, are given in Table 6.1. The relative error between our•

9Table I. Average PDg Performance for Rigid Models with Two Different
Free-configuration Selection Methods.
Benchmark

Motion Coherence
Timing# of Iterations

Centroid Difference
Timing# of Iterations

Bunny/Dragon
Spoon/Cup

22.8 msec
7.76 msec

59.2 msec
24.1 msec

2.14
2.35

6.5
7.3

Table II. Performance Comparison between Our Method and[Zhang et al. 2007a].
Benchmark

Ours

Zhang et al.

Relative Error

Pawn/CAD
Hammer/Bumpy-Sphere
L-shape/Torusknot

11.1 msec
4.9 msec
3.2 msec

80.2 msec
57.7 msec
186 msec

0.21
0.02
0.32

and Zhang et al.’s method was computed as:e≡|PD1g − PD2g |max PD1g , PD2g

,

(31)

where, PD1g and PD2g are the results of our and Zhang et al.’s,
respectively. In the experiments, our algorithm outperforms their
method by a factor of 7 ∼ 60 with small average relative-error differences; even in the case of the L-shape/Torusknot, our results
generate smaller object norm values than [Zhang et al. 2007a] formore than 65% of simulation frames.

Fig. 9. Rigid Benchmarking Models with Triangle Count. CAD (2.4K),
Spoon (1.34K), Dragon (174K), Bunny (40K), Torusknot (3K), L-shape(20), Cup (1K), Pawn (0.3K), Bumpy-Sphere (2.9K), Hammer (1.7K).

PDg Computation for Articulated Models We use three testcases to benchmark the performance of our algorithm, also shownin Fig.14.
(1) Puma1 (Fig. 14 (a)): a Puma robot with 8 links, 0.9K triangles
and 9-DoF translates along a Beetle car, consisting of 4K triangles, and collides with it (but no self-collision for the robot).
(2) Puma2 (Fig. 14 (b)): the same Puma robot rotates its upperarm to create self-collision with itself, and also collides withthe car.
(3) Hand (Fig. 14 (c)): a hand model (25 links, 6K triangles, and25-DoF) translates from top to down, while colliding with abottle (0.7K triangles).
The average statistics of our algorithm are summarized in Table.III.

6.2

Quantitative Error Analysis

To the best of our knowledge, there exists no algorithm to exactly compute generalized penetration depth (PDg ) for either generalrigid polygonal or articulated models, due to the high dimensionality of the contact surfaces. Thus, designing and performing comparative evaluations for general polygonal models are nearly impossible. However, in some limited cases, one can design a samplingbased, brute force method to rather accurately approximate PDg for

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

10•Tang and Kim(a) Spoon/Cup(b) Bunny/Dragon(c) Placement(d) Robot under the Bridge(e) Grasping

Fig. 10. Interactive Generalized Penetration Depth (PDg ) Computation (a, b) As a result of PDg computation for rigid models, the red, colliding objects
are rigidly transformed to the green ones, just in contact with obstacles. PDg computation for articulated models and application to motion planning (c, d)and grasping (e). Our algorithm can compute PDg for these challenging benchmarks at interactive rates. The models in (c), (d) and (e) were obtained fromthe GAMMA and KIT research groups, respectively, under permission.

Table III. Average PDg Performance for
Articulated Models.
BenchmarkTimingNumber of Iterations

Puma1
Puma2
Hand

27.9 msec
14.9 msec
71.6 msec

3.1
2.3
1.9

long each dimension. In this experiment, the average relative error
between the brute force method and ours using Eq. 31 is 0.0165.

a rigid body. We compare the results of our PDg algorithm against
this brute force algorithm. For an articulated body, we compare ourmethod with an optimal result that can be analytically calculated
for a special case.
Rigid body As shown in Fig. 11, we compare the results of our algorithm against the brute force method for a relatively simple case:a moving, polygonal spoon model (displayed in red color) penetrates into a stationary box model, and the PDg between them iscalculated repeatedly by translating the spoon model (displayed incyan color). In the figure, the red, green and blue arrows represent
the X, Y and Z axes of the local frame attached to the spoon model. For the spoon model, Izz is ten times greater than Iyy . Thus, we
can just consider only rotation with respect to the Y axis to get the
optimum of PDg . Further, since any change in the rotation angle
about X axis will induce more rotation about the Y axis, to get acontact configuration (i.e. the cyan configuration), it will result ina higher value of object norm; thus we need not consider rotation
with respect to the X axis. Then, the total DoF of rotation that weneed to consider can be decreased from three to one.
Furthermore, since the box has a planar surface (and thus a planar
contact), the translational component of PDg result should be on
the X-Z plane defined by the origin of the penetrated configurationand the rotational axis (Y axis in the case). The reason is that, for
any contact configuration whose origin is not on the X-Z plane, one
can always project it back to the plane with a smaller object norm.
Then, the DoF of translational component in PDg can be reduced
from three to two. Also, note that this argument is irrelevant to theorientation of the spoon model.
Thus, the final DoF that we need to consider to get an optimal valueof PDg is three (one rotation + two translations). In this case, one
can rather accurately approximate the PDg by sampling the contact
space as follows: we get the translational components of PDg by
sampling the X-Z plane of the red spoon and get the contact configurations by sampling the rotational angle about the Y axis. The
contact configuration with a minimum object norm is chosen as the
PDg configuration. In our implementation, we sample 300 times a-

Fig. 11. Sampling-based Method to Accurately Approximate the Optimal PDg . A spoon model (red) penetrates into a stationary box, and the
PDg between them (cyan) is calculated repeatedly, by placing the spoon
model at different configurations.

Articulated body We analytically calculate the optimal result for a
special case: a planar articulated model composed of two stick-likelinks with revolute joint penetrates an infinitely large ground obstacle, as illustrated in Fig.12. Here, one end of the root link is fixedand displaced from the ground by d, and the joint parameters are[θ1 , θ2 ], the length of each link is l. In the experiment, we change d
to generate different generalized penetration depth values.

d

1 l

l

2
Obstacle

Fig. 12. A Special Case to Calculate the Optimal PDg . The blue articulated model with two links is placed at a configuration {θ1 , θ2 }, displaced
from the red one at a rest configuration.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

•Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Norm

11

In this simple case, the object norm of the model can be analyticallycalculated as  3
  3
1
θ1 8l
θ2 l
σ=
4 sin2+ 4 sin2
(32)
V
2
3
2 3
where l cos θ1 + l cos (θ1 + θ2 ) = d.
By differentiating Eq. 32, the optimal PDg value can be obtained.
In Fig.13, we show the comparisons between the optimal and our
PDg results as d changes from 1.0l to 1.9l . The mean squared error
is 0.0027, and our result is very close to the optimum.
0.7

Our method
Optimum

0.6

PDg value

0.5
0.4
0.3
0.2
0.1
0
1

1.1

1.2

1.3

1.4

d

1.5

1.6

1.7

1.8

1.9

(b) Puma2(c) HandFig. 14. PDg Computation for Articulated Models. The models in (a) and(b) are included as part of the open-sourced, motion planning kit (MPK).
The models in (c) were obtained from the KIT research group under permission.

algorithm (BiRRT-PDg ) are shown in Table IV. These results wereobtained by running each of the three planners ten times, and averaging their timings, and relevant sampling parameters were kept thesame during the experiments. The new planner using our PDg algorithm outperforms all other algorithms by a factor of 1.8, because
the retracted samples generated by our algorithm are closer to thein-collision samples than others. Moreover, the result of planned
paths using our algorithm look visually more plausible than others (as shown in the accompanying video), since our distance metricbased on object norm discourages the large motion of an object that
requires more energy.
Table IV. Average Performance of Different
Motion Planners.

Fig. 13. Comparisons between the Optimum and Our PDg Method.

6.3(a) Puma1

Applications

In this section, we demonstrate the use of our PDg algorithm inthree different applications.
Motion Planning In sampling-based motion planning, the narrowpassage problem is a known hard problem, and some efficient algorithms based on retraction have been proposed in the literature[Redon and Lin 2005; Zhang and Manocha 2008; Pan et al. 2010].
During the sampling process, these algorithms attempt to retract anin-collision sample toward a free configuration, yielding a contact
sample close to the in-collision sample [Redon and Lin 2005; Panet al. 2010]. We use our PDg algorithm to generate contact samples
for retraction-based planners.
More specifically, we integrate our PDg algorithm into BiRRT[Kuffner and LaValle 2000], by retracting the in-collision sample
q to the contact space, using our PDg algorithm; we use a free nodein the RRT tree nearest to q as an initial, free configuration q f for
PDg optimization. To measure the performance of the BiRRT algorithm with PDg (BiRRT-PDg in Table IV), we use two planningscenarios, similarly used in [Pan et al. 2010]:(1) Placement (Fig.10(c)): A human-like figure with 15 links and
41 DoFs picks a tool from a desk and places it on a book shelf.
The environment including the figure consists of 52K triangles.
(2) Bridge (Fig.10(d)): A robot with 13 links and 40 DoFs moves
through a narrow hole under the bridge.The environment including the robot consists of 31K triangles.
In Figs .10 (c) & (d), the red and green figures/robots represent
their start and goal configurations, respectively. The performancesof different retraction-based algorithms using BiRRT [Kuffner and
LaValle 2000], BiRRT-R [Pan et al. 2010], and BiRRT using our

Benchmark

BiRRT

BiRRT-R

BiRRT-PDg

Placement
Bridge172.8s78.3s154.2s

42.1s
85.3s

Physically-based Animation There can be many ways to utilize
PDg algorithms in physically-based animation; e.g. [Ortega et al.
2007; Kim et al. 2003]. Here, we propose a novel way to use PDg
for rigid or articulated-body dynamics. In these applications, as an
initial condition for simulation, all participating objects should bein a collision-free state. However, doing so is non-trivial as well as
manually intensive for a user to enforce such non-penetration constraints, especially for a complicated articulated character, since it
requires repetitive invocation of collision detection until a satisfactory, collision-free pose is obtained for the character. In this case,
we propose to use our PDg algorithm; the user will sketch a desired
but possibly colliding pose for an articulated character and then call PDg to fix the pose so that it is collision-free as well as close to
the intended pose. In Fig.15, we created two possible scenarios forarticulated-body dynamics: a figure sitting on a bench and sitting
inside a cup. In the figure, a user create a desired pose in red, but incollision, and the PDg algorithm can generate a collision-free posein green. The bench and cup benchmark takes 85.2 msec and 87.2
msec to compute PDg . A known physics engine such as Bullet2
does not work when a character is given as the red pose, but withthe green poses, physics simulations work.
Grasping In data-driven grasp planning, it is important to quickly
find valid contacts for a novel object from the grasping pose of asimilar object from a database [Goldfeder et al. 2009]. Moreover,since the grasping pose in the database is optimal in terms of grasping quality, the contact pose found for the novel object should beclose to the optimal one if it is collision-free.
2 http://bulletphysics.org

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

12•Tang and Kim—We locally linearize the nonlinear contact surface to turn the nonlinear optimization problem of PDg calculation to a linear one.
This simple linearization works well in practice, since the interpenetration between objects is mostly local and shallow in mostapplications.
—Our distance metric (or objective function), the object norm, hasa closed form of a polynomial function, which makes the optimization rather easy.
(a) Bench(b) CupFig. 15. Finding Collision-free Configurations for Articulated-body Dynamics.

For such a purpose, we can use our PDg algorithm. For instance, inFig. 10 (e), to grasp the red cup model (middle), a similar green cup
model and its optimal pose are searched from the database (right).
Then, the pose from the database is used as a pre-grasp for the redcup model (middle), but this pose will generate collision. Then, ourPDg algorithm is applied to generate a valid pose (left), which is
close to the optimal pose for the red cup. Our PDg algorithm takes
148 msec to compute a valid, collision-free pose for the hand/cup
model.

6.4

Discussions

In this section, we show the computational complexity of our PDg
algorithm, compare our algorithm against other known algorithms,
and explain some limitations of our algorithm.
Computational Complexity In the case of PDg for a rigid model,
its computational cost is similar to [Je et al. 2012], except for someconstant difference to set up the LLCS instead of LCS. So, we refer
to [Je et al. 2012] for more detailed analysis.
The main computational cost dominating our PDg algorithm foran articulated model includes free-configuration selection, contactspace projection and constrained optimization. Except for successive perturbation, our free-configuration selection scheme is O(n),where n is the number of links in A . In case of successive perturbation, each invocation of perturbation requires repeated collisionchecks for a small range of parameter intervals. Thus, the complexity of successive perturbation is TSP = O(pc nTCD ), where pc
is the maximum number of collision checks for a parameter interval, defined by the user, and TCD is the cost for collision detection for articulated models [Schmidl et al. 2004]. The contactspace projection requires CCD operations n times in the worst
case, since kLk ≤ n and kTk ≤ n in Alg.2. Thus, the complexity of CSP is TCSP = O(nTCCD ) where TCCD is the cost to perform
CCD like [Zhang et al. 2007b]. The constrained optimization isdominated by the LCP solver. Since we use the Gauss-Seidel (GS) method to solve an LCP, its complexity is TLCP = O(mGS c2 )where mGS is the number of G-S iterations and c is the maximum number of contact points per each G-S iteration. Therefore,
the complexity of our PDg algorithm for an articulated model is
T = max (O(n), TSP ) + N(TCSP + TLCP ), where N is the number of
iterations for constrained optimization, which is typically 2∼3 inpractice.
Despite the worst-case computational complexity of our algorithm,
our algorithm performs at interactive rates in practice, due to thefollowing reasons:—Our strategies for free-configuration selection work well in practice, in particular when the motion coherence is present in a simulated environment, which is often the case for physically-based
animation and sampling-based motion planning.
Convergence Analysis The convergence of our algorithm is guaranteed, due to the following reasons:—The constrained optimization in our case was formulated as an
LCP, and the PGS was employed to solve the LCP; refer to Eq.
12. In this case, the PGS is convergent, since JM−1 JT in Eq. 12
is an SPD matrix. So is for articulated model, since C0 H−1 C0 T
is SPD in Eq. 30.
—The contact space projection performs pretty much like CCD algorithms based on conservative advancement. Its convergencecan be supported using the same argument presented in earlier
CCD work, like [Zhang et al. 2007b].
—At every iteration, the constrained optimization always attempts
to decrease the object norm value, and so does the contact space
projection. In particular, the constrained optimization step monotonically decreases the object norm at every iteration; otherwise,
we stop optimizing it (i.e. we found a local optimum).
Our algorithm terminates after a finite number of steps, since thecomplexity of contact space is bounded above.
Comparisons We compare the main differences of our algorithmwith [Zhang et al. 2007a; Zhang et al. 2007a; Nawratil et al. 2009;
Je et al. 2012] in terms of PDg computation for rigid models.
In [Zhang et al. 2007a], the PDg problem is posed as a containment
problem, and linear programming is employed to solve it. Theyuse a distance metric called Dg defined using trajectory length,
and their approach is limited to a certain class of objects. Moreover, their method can compute only bounds of PDg . Zhang et al.’s
method [Zhang et al. 2007a] also relies on the object norm andconstrained optimization to find a solution for PDg . However, themain difference between our method and [Zhang et al. 2007a] is the
optimization step. Zhang et al.’s algorithm relies on random sampling on contact space to optimize the object norm, whereas our
method uses an LCP-based formulation that does not need sampling. This is the main reason why our method outperforms theirs
by an order of magnitude. Our algorithm shares some similarity
with [Nawratil et al. 2009] in that both of the methods are based on
iterative, constrained optimization. We can also prove that the metric S in [Nawratil et al. 2009] is equivalent to object norm undersome condition, as shown in Appendix B. However, this algorithm is relatively slow in terms of performance, since it relies on an
adjustment step based on distance fields and free-configuration selection scheme that are costly, and thus it may not be suitable forinteractive application. Finally, [Je et al. 2012] also uses contactspace projection to iteratively optimization the PD. However, their
method is applicable only to linear contact-space, and thus translational penetration depth, and the underlying distance metric is
Euclidean.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

•Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Norm

It is not clear whether any of the above work is applicable to articulated models.
Limitations Our algorithm has some limitations. Our algorithm
has no guarantee that the found PDg solution is globally optimal,and our algorithm can be stuck in local minima. Moreover, the result of our algorithm highly depends on the selection of an initialcollision-free configuration. However, in a practical setting such as
our benchmarking examples, our algorithm generates a fairly goodPDg solution, so that many applications can benefit from the result.
Our algorithm may require a large number of iterations to converge
for complicated, rigid or articulated models with high DoFs and alarge number of contacts. Currently, our PDg algorithm cannot handle an articulated model with a closed loop, though such a model is
not popularly used in graphics.

A.

CONTACT-SPACE REALIZATION

Theorem 1 For given, colliding articulated models A and B
with no closed loops and A being initially placed at an incollision configuration qin , if q is a configuration of A that realizes
PDg (A , B), then q is on the contact space induced by A and B.
PROOF We prove it by contradiction. Suppose that q realizing PDg (A , B) is not on the contact space; then, q must be infree space. Let us change the parameters for each joint sequentially, and first consider the base joint J0 of A . In this case, wecreate a new configuration q0 for A , where q0 (J0 ) = q(J0 ) andq0 (Ji ) = q0 (Ji ), i = 1, · · · , n − 1, as illustrated in Fig.16(a). By thedefinition of object norm for a single joint as Eq. 13, we have
σ0 (qin , q0 ) = σ0 (qin , q).

7.

13

(33)

Then, q0 can have two cases:

CONCLUSIONSIn this paper, we present a real-time algorithm to compute generalized PD for rigid and articulated models. The generalized PD
is formulated as a minimum distance query on the contact space using the object norm. Our algorithm relies on iterative constrained optimization on the contact space consisting of four steps:
free-configuration selection, CSP, constrained optimization and reprojection. We have also applied our algorithm to different applications including motion planning, physically-based animation andgrasping.
Future Work There are a few interesting future directions for ourwork. Extension of the generalized PD computation to multiple
bodies or deformable bodies, including its proper definition will
be interesting. Many other applications can potentially get benefit
from the generalized PD. For instance, the generalized PD (in particular, the angular term) may be integrated into constraint-based
dynamics, for a better convergence for a contact-constraint solver,
especially when the rotational motion in simulation is large. The
angular term may be effectively utilized in 6-DOF haptics, to generate proper torque feedback. GPU-based or parallel implementations of our algorithm will be also an interesting research direction
to pursue. Currently, our algorithm can calculate PDg for a pair ofmoving and stationary objects. It will be interesting to find PDg for
N moving, overlapping bodies to separate them simultaneously.

(1) q0 in free space. In this case, A can rigidly move from qin
to q0 , since qin and q0 are different only for the base joint J0 .
As proved by [Zhang et al. 2007a], there must be a contact configuration qc for A when rigidly moving from qin to q0 ,which makes σ0 (qin , qc ) < σ0 (qin , q0 ). According to the definition of object norm for the whole articulated model (Eq.19)
and according to Eq.33, and because σi (qin , qc ) = 0, for i =1, · · · , n − 1, we have σ (qin , qc ) = σ0 (qin , qc ) < σ0 (qin , q0 ) =σ0 (qin , q) < σ (qin , q), which contradicts to our assumptionthat q realizes PDg (A , B).
(2) q0 is in obstacle space. In this case, we modify q0 suchthat q0 (J0 ) = q(J0 ), q0 (J1 ) = q(J1 ), and q0 (Ji ) = qin (Ji ), i =2, · · · , n, as illustrated in Fig.16(b). Similarly to the previous
case, the new q0 can be either in free space or in obstacle space. If q0 is in free space, we can find a contact configuration
qc such that σ1 (qin , qc ) < σ1 (qin , q0 ), and σi (qin , qc ) = 0, i =2, · · · , n − 1. Then, σ (qin , qc ) < σ (qin , q), which is a contradiction. If q0 is again in obstacle space, a recursive argument can be made as follows. Since q is in free space, we must
be able to find a joint Jk and a free configuration q0 such that
q0 (Ji ) = q(Ji ), i ≤ k and q0 (J j ) = qin (J j ), j > k. Then, a contact
configuration qc can be found so that σk (qin , qc ) < σk (qin , q0 ),which makes σ (qin , qc ) < σ (qin , q). q

ACKNOWLEDGEMENTS

q

q'

q'

This research was supported in part by NRF in Korea(No.2012R1A2A2A01046246, No.2012R1A2A2A06047007). We
thank Liangjun Zhang for helpful initial discussions about thisproject.

qc

qc

J0

qin
obstacleJ0

qin
obstacle

APPENDIX(a)

We prove the contact-space realization theorem for articulated
models, and also show the relationship between object norm [Kazerounian and Rastegar 1992] and metric S [Nawratil et al. 2009] fora rigid model.

(b)Fig. 16. The intermediate configuration q0 . A rigid motion from q to q0
can be created, since the joint parameters are the same except for the base.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

14

B.

•

Tang and Kim

EQUIVALENCE RELATIONSHIP BETWEEN METRIC

By comparing Eq. 38 with Eq. 39, we can see that

S AND OBJECT NORM
Give a rigid model A , the metric S is defined as the squared distance between six special points of an object A at two different
configurations q0 , q1 [Nawratil et al. 2009]:

1 2
S = σP (q0 , q1 ) , when N = 6.
6

(40)

Thus, when N = 6, S2 is equivalent to the discrete the object normof the point set, defined by the vertices of A . In this sense, themetric S can be considered a discrete version of the object norm.

6

S (q0 , q1 )2 = ∑ ksi (q0 ) − si (q1 )k2 ,(34)
REFERENCES

i=1

where the six special points si , i = 1, ...6 are defined as:qq
si = bx + λ2i di , si+3 = bx − λ2i di , i = 1, 2, 3.

(35)

In Eq. 35, bx is the center of mass of A , λi and di are the eigenvalue
and the corresponding unit eigenvector of a covariance matrix Dx =NT

∑ (xi − bx ) (xi − bx ) , where xi , i = 1, ...N are the vertices of ai=1

triangle mesh of A .

If we define the world frame at the center of mass and principle axesof A , we can relate the metric S to the object norm, by expressing
the metric S in terms of q, which is used in Eq.1 to define the object
norm for a rigid model: We set the affine map for a configuration
qA ,1 as A = [Ta , Ra ], where Ta is a translation vector, and Ra =[a1 , a2 , a3 ] is a rotation matrix. According to [Hofer and Pottmann2004], the metric S also can be expressed as:
3

S2 = 6c20 + 2 ∑ fi2 c2i ,i=1where c0 = Ta = [q4 , q5 , q6 ]T , fi , i = 1, 2, 3 are the eigenvaluesof Dx , ci , i = 1, 2, 3 are the vectors of Ra − E. We now use q =[q1 , q2 , q3 , q4 , q5 , q6 ] to represent metric S. Since
T 
 
1
1
c21 = a1 −  0  a1 −  0 
00
= (a11 − 1)2 + a212 + a213
= 2 − 2a11

= 2 1 − q20 − q21 + q22 + q23
22= 4 q2 + q3 ,where ai j is an element of Ra , and because


N
∑ x2
0 + I0 − I0
Iyy
zz
xx
=
.
f12 = i=1
2
4

(36)

(37)

Note that the inertia tensor I 0 used in Eq.37 is computed from only
the vertices of A , not from the entire body of A .
We can derive similar equations like Eq.36 and Eq.37 for ci , i = 2, 3
and fi , i = 2, 3. Then, the metric S can be represented as:




0 2
0 2
0 2
S2 = 4 Ixx
q1 + Iyy
q2 + Izz
q3 + 6 q24 + q25 + q26 .
(38)
Now, let us consider a discrete version of the object norm for apoint set P = {pi , i = 1, . . . , N} consisting of A ’s vertices. Then,σP (q0 , q1 )= N1 ∑ (x (q0 ) − x (q1 ))2
x∈P
0 q2 + I 0 q2 + I 0 q2 + q2 + q2 + q2 .
= N4 Ixx
yy 2
zz 3
1
4
5
6

(39)

P. K. Agarwal, L. J. Guibas, S. Har-peled, A. Rabinovitch, and M. Sharir.
2000. Penetration depth of two convex polytopes in 3D. Nordic J. Computing 7 (2000), 227–240.
J. Allard, F. Faure, H. Courtecuisse, F. Falipou, C. Duriez, and P. G. Kry.
2010. Volume contact constraints at arbitrary resolution. ACM Trans.
Graph. 29, Article 82 (July 2010), 10 pages. Issue 4.
J. Baumgarte. 1972. Stabilization of constraints and integrals of motion
in dynamical systems. Computer Methods in Applied Mechanics andEngineering 1, 1 (1972), 1 – 16.
G. Bergen. 2001. Proximity Queries and Penetration Depth Computation
on 3D Game Objects. Game Developers Conference (2001).
S. Cameron. 1997. Enhancing GJK: Computing Minimum and Penetration Distance between Convex Polyhedra. Proceedings of International
Conference on Robotics and Automation (1997), 3112–3117.
S. A. Cameron and R. K. Culley. 1986. Determining the minimum translational distance between two convex polyhedra. In Proc. of IEEE Inter.
Conf. on Robotics and Automation. 591–596.
E. Catto. 2005. Iterative Dynamics with Temporal Coherence. In Game
Developer Conference. San Francisco, USA.
D. Dobkin, J. Hershberger, D. Kirkpatrick, and S. Suri. 1993. Computing
the intersection-depth of polyhedra. Algorithmica 9 (1993), 518–533.
S. Fisher and M. C. Lin. 2001. Fast Penetration Depth Estimation for Elastic
Bodies Using Deformed Distance Fields. Proc. of IEEE/RSJ International Conference on Intelligent Robots and Systems 2001 (2001).
C. Goldfeder, M. Ciocarlie, H. Dang, and P. K. Allen. 2009. The Columbia
Grasp Database. In International Conference on Robotics and Automation.
P. Hachenberger. 2009. Exact Minkowksi Sums of Polyhedra and Exact andEfficient Decomposition of Polyhedra into Convex Pieces. Algorithmica
55, 2 (2009), 329–345.
G. D. Hart, B. M. Riviere, and G. D. Hart. 2002. A constraint-stabilized
time-stepping approach for rigid multibody dynamics with joints, contact
and friction. In International Journal for Numerical Methods in Engineering. 2335–2371.
M. Hofer and H. Pottmann. 2004. Energy-Minimizing Splines in Manifolds.
Proc. of ACM SIGGRAPH (2004).
K. Hoff, A. Zaferakis, M. Lin, and D. Manocha. 2002. Fast 3d Geometric
Proximity Queries Between Rigid and Deformable Models Using Graphics Hardware Acceleration. Technical Report TR02-004. Department ofComputer Science, University of North Carolina.
C. Je, M. Tang, Y. Lee, M. Lee, and Y. Kim. 2012. PolyDepth: Real-time
Penetration Depth Computation using Iterative Contact-Space Projection.
ACM Trans. Graph. 31, Article 5 (January 2012). Issue 1.
F. Jourdan, P. Alart, and M. Jean. 1998. A Gauss-Seidel like algorithm
to solve frictional contact problems. Computer Methods in Applied Mechanics and Engineering 155, 1-2 (1998), 31–47.
K. Kazerounian and J. Rastegar. 1992. Object norms: A class of coordinate
and metric independent norms for displacement, In in Flexible Mechanism, Dynamics and Analysis: ASME Design Technical Conference. 22nd
Biennial Mechanisms Conference, G. K. et al, Ed. 47, 1-2 (1992), 271–
275.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Interactive Generalized Penetration Depth Computation for Rigid and Articulated Models using Object Norm
Y. J. Kim, M. C. Lin, and D. Manocha. 2004. Incremental Penetration
Depth Estimation between Convex Polytopes using Dual-Space Expansion. IEEE Trans. on Visualization and Computer Graphics 10, 1 (2004),
152–164.
Y. J. Kim, M. A. Otaduy, M. C. Lin, and D. Manocha. 2002. Fast Penetration Depth Computation for Physically-based Animation. Proc. of ACM
Symposium on Computer Animation (2002).
Y. J. Kim, M. A. Otaduy, M. C. Lin, and D. Manocha. 2003. Six-Degreeof-Freedom Haptic Rendering Using Incremental and Localized Computations. Presence 12, 3 (2003), 277–295.
J. J. Kuffner and S. M. LaValle. 2000. RRT-Connect: An efficient approach
to single-query path planning. In Proc. IEEE International Conferenceon Robotics and Automation.
E. Larsen, S. Gottschalk, M. C. Lin, and D. Manocha. 2000. Fast Proximity
Queries with Swept Sphere Volumes. IEEE International Conference on
Robotics and Automation (ICRA) (2000), 3719–3726.
J. Lien. 2008. Covering Minkowski Sum Boundary Using Points with Applications. Computer Aided Geometric Design 25, 8 (November 2008),
652–666.
J. Lien. 2009. A Simple Method for Computing Minkowski Sum Boundary
in 3D Using Collision Detection. In Algorithmic Foundation of Robotics VIII, G. Chirikjian, H. Choset, M. Morales, and T. Murphey (Eds.).
Springer Tracts in Advanced Robotics, Vol. 57. Springer Berlin / Heidelberg, 401–415.
B. Mirtich. 1996. Impulse-based Dynamic Simulation of Rigid Body Systems. Ph.D. Dissertation. University of California, Berkeley.
M. Moore and J. Wilhelms. 1988. Collision Detection and Response for
Computer Animation. In Computer Graphics (SIGGRAPH ’88 Proceedings), John Dill (Ed.), Vol. 22. 289–298.
G. Nawratil, H. Pottmann, and B. Ravani. 2009. Generalized penetration
depth computation based on kinematical geometry. Comput. Aided Geom. Des. 26 (May 2009), 425–443. Issue 4.
C.J. Ong. 1993. Penetration distances and their applications to path planning. Ph.D. Dissertation. Michigan Univ., Ann Arbor.
C. J. Ong and E.G. Gilbert. 1996. Growth Distances: New Measures for
Object Separation and Penetration. IEEE Transactions on Robotics andAutomation 12, 6 (1996), 888–903.
M. Ortega, S. Redon, and S. Coquillart. 2007. A Six Degree-of-Freedom
God-Object Method for Haptic Display of Rigid Bodies with Surface
Properties. IEEE Transactions on Visualization and Computer Graphics 13 (May 2007), 458–469. Issue 3.
M. A. Otaduy, R. Tamstorf, D. Steinemann, and M. Gross. 2009. Implicit
Contact Handling for Deformable Objects. Computer Graphics Forum(Proc. of Eurographics) 28, 2 (Apr 2009).
J. Pan, L. Zhang, and D. Manocha. 2010. Retraction-Based RRT Planner
for Articulated Models. In IEEE International Conference on Robotics
and Automation (ICRA).
S. Redon, A. Kheddar, and S. Coquillart. 2002. Gauss’ least constraints
principle and rigid body simulations. In In proceedings of IEEE International Conference on Robotics and Automation. 11–15.
S. Redon and M. Lin. 2005. Practical local planning in the contact space.
In IEEE International Conference on Robotics and Automation(ICRA).
S. Redon and M. C. Lin. 2006. A Fast Method for Local Penetration Depth
Computation. Journal of graphics tools 11, 2 (2006), 37–50.
D. C. Ruspini and O. Khatib. 1997. Collision/Contact Models for the Dynamic Simulation of Complex Environments. In In 9th International Symposium of Robotics Research. Snowbird, USA.
H. Schmidl, N. Walker, and M. C. Lin. 2004. CAB: Fast update of OBB trees for collision detection between articulated bodies. Journal ofGraphics Tools 9 (2004), 1–9.

•

15

A. Sud, N. Govindaraju, R. Gayle, I. Kabul, and D. Manocha. 2006. Fast
proximity computation among deformable models using discrete Voronoi
diagrams. In ACM SIGGRAPH. ACM New York, NY, USA, 1144–1153.
M. Tang, Y. J. Kim, and D. Manocha. 2009a. C2 A: controlled conservativeadvancement for continuous collision detection of polygonal models. In
IEEE International Conference on Robotics and Automation. 356–361.
M. Tang, M. Lee, and Y. J. Kim. 2009b. Interactive Hausdorff distance computation for general polygonal models. ACM Trans. Graph. 28, Article
74 (July 2009), 9 pages. Issue 3.
B. Wang, F. Faure, and D. K. Pai. 2012. Adaptive Image-based Intersection
Volume. ACM Trans. Graph. (Proc. SIGGRAPH) 31, 4 (2012).
R. Weller and G. Zachmann. 2009. Inner sphere trees for proximity andpenetration queries. In Proceedings of Robotics: Science and Systems.
Seattle, USA.
L. Zhang, Y. Kim, and D. Manocha. 2007a. A fast and practical algorithm for generalized penetration depth computation. In Proceedings of
Robotics: Science and Systems. Atlanta, USA.
L. Zhang, Y. J. Kim, and D. Manocha. 2007b. C-DIST: Efficient DistanceComputation for Rigid and Articulated Models in Configuration Space. In
ACM Solid and Physical Modeling Conference (SPM07). Beijing, China.
L. Zhang, Y. J. Kim, and D. Manocha. 2008. Efficient Cell Labelling andPath Non-existence Computation using C-obstacle Query. The International Journal of Robotics Research 27, 11-12 (Nov-Dec 2008), 1246–
1257.
L. Zhang, Y. J. Kim, G. Varadhan, and D. Manocha. 2007a. Generalized
penetration depth computation. Computer-Aided Design 39, 8 (2007),
625–638.
L. Zhang and D. Manocha. 2008. An Efficient Retraction-based RRT Planner. In IEEE International Conference on Robotics and Automation (ICRA). 3743–3750.
X. Zhang, S. Redon, M. Lee, and Y. J. Kim. 2007b. Continuous collision detection for articulated models using Taylor models and temporal culling.
ACM Trans. Graph. 26, Article 15 (July 2007). Issue 3.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.