Skip to content

Commit c3c47c7

Browse files
authored
[FEATURE] Implemented GJK collision detection algorithm (#1213) (#1213)
1 parent 01ab1ae commit c3c47c7

File tree

12 files changed

+3428
-145
lines changed

12 files changed

+3428
-145
lines changed

genesis/assets/xml/tet_ball.xml

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
<mujoco model="tet_ball">
2+
<compiler angle="degree"/>
3+
<option gravity="0 0 -9.8100000000000005"/>
4+
<default>
5+
<default class="/">
6+
<joint damping="0.001"/>
7+
<geom friction="1 0.0050000000000000001 0.0001" density="500"/>
8+
</default>
9+
</default>
10+
11+
<asset>
12+
<texture builtin="gradient" height="100" rgb1="1 1 1" rgb2="0 0 0" type="skybox" width="100"/>
13+
<texture builtin="flat" height="1278" mark="cross" markrgb="1 1 1" name="texgeom" random="0.01" rgb1="0.8 0.6 0.4" rgb2="0.8 0.6 0.4" type="cube" width="127"/>
14+
<texture builtin="checker" height="100" name="texplane" rgb1="0 0 0" rgb2="0.8 0.8 0.8" type="2d" width="100"/>
15+
<material name="MatPlane" reflectance="0.5" shininess="1" specular="1" texrepeat="60 60" texture="texplane"/>
16+
<mesh name="tet" file="tet.obj"/>
17+
</asset>
18+
19+
<worldbody>
20+
<geom name="tet1" contype="1" conaffinity="1" pos="0 0 0" type="mesh" rgba="0 0 1 1" mesh="tet"/>
21+
22+
<body name="sphere1" pos="0 0 1.2">
23+
<joint name="root1" class="/" type="free"/>
24+
<geom name="seg1_geom" class="/" type="sphere" size="0.1"/>
25+
</body>
26+
27+
<body name="sphere2" pos="0.3 0 1.2">
28+
<joint name="root2" class="/" type="free"/>
29+
<geom name="seg2_geom" class="/" type="sphere" size="0.1"/>
30+
</body>
31+
32+
<body name="sphere3" pos="0.3 0.29 1.2">
33+
<joint name="root3" class="/" type="free"/>
34+
<geom name="seg3_geom" class="/" type="sphere" size="0.1"/>
35+
</body>
36+
</worldbody>
37+
</mujoco>

genesis/assets/xml/tet_capsule.xml

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
<mujoco model="tet_ball">
2+
<compiler angle="degree"/>
3+
<option gravity="0 0 -9.8100000000000005"/>
4+
<default>
5+
<default class="/">
6+
<joint damping="0.001"/>
7+
<geom friction="1 0.0050000000000000001 0.0001" density="500"/>
8+
</default>
9+
</default>
10+
11+
<asset>
12+
<texture builtin="gradient" height="100" rgb1="1 1 1" rgb2="0 0 0" type="skybox" width="100"/>
13+
<texture builtin="flat" height="1278" mark="cross" markrgb="1 1 1" name="texgeom" random="0.01" rgb1="0.8 0.6 0.4" rgb2="0.8 0.6 0.4" type="cube" width="127"/>
14+
<texture builtin="checker" height="100" name="texplane" rgb1="0 0 0" rgb2="0.8 0.8 0.8" type="2d" width="100"/>
15+
<material name="MatPlane" reflectance="0.5" shininess="1" specular="1" texrepeat="60 60" texture="texplane"/>
16+
<mesh name="tet" file="tet.obj"/>
17+
</asset>
18+
19+
<worldbody>
20+
<geom name="tet1" contype="1" conaffinity="1" pos="0 0 0" type="mesh" rgba="0 0 1 1" mesh="tet"/>
21+
22+
<body name="capsule1" pos="0 0 1.2">
23+
<joint name="root1" class="/" type="free"/>
24+
<geom name="seg1_geom" class="/" type="capsule" size="0.02 0.05"/>
25+
</body>
26+
27+
<body name="capsule2" pos="0.25 0 1.2">
28+
<joint name="root2" class="/" type="free"/>
29+
<geom name="seg2_geom" class="/" type="capsule" size="0.04 0.1"/>
30+
</body>
31+
32+
<body name="capsule3" pos="0.5 0 1.2">
33+
<joint name="root3" class="/" type="free"/>
34+
<geom name="seg3_geom" class="/" type="capsule" size="0.1 0.05"/>
35+
</body>
36+
</worldbody>
37+
</mujoco>

genesis/assets/xml/tet_tet.xml

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
<mujoco model="two_tet">
2+
<compiler angle="degree"/>
3+
<option gravity="0 0 -9.8100000000000005"/>
4+
<default>
5+
<default class="/">
6+
<joint damping="0.001"/>
7+
<geom friction="1 0.0050000000000000001 0.0001" density="500"/>
8+
</default>
9+
</default>
10+
11+
<asset>
12+
<texture builtin="gradient" height="100" rgb1="1 1 1" rgb2="0 0 0" type="skybox" width="100"/>
13+
<texture builtin="flat" height="1278" mark="cross" markrgb="1 1 1" name="texgeom" random="0.01" rgb1="0.8 0.6 0.4" rgb2="0.8 0.6 0.4" type="cube" width="127"/>
14+
<texture builtin="checker" height="100" name="texplane" rgb1="0 0 0" rgb2="0.8 0.8 0.8" type="2d" width="100"/>
15+
<material name="MatPlane" reflectance="0.5" shininess="1" specular="1" texrepeat="60 60" texture="texplane"/>
16+
<mesh name="tet" file="tet.obj"/>
17+
<mesh name="small_tet" file="tet.obj" scale="0.1 0.1 0.1"/>
18+
</asset>
19+
20+
<worldbody>
21+
<geom name="tet1" contype="1" conaffinity="1" pos="0 0 0" type="mesh" rgba="0 0 1 1" mesh="tet"/>
22+
23+
<body name="tet2" pos="0 0 1.5">
24+
<joint name="root1" class="/" type="free"/>
25+
<geom name="tet2" contype="1" conaffinity="1" pos="0 0 0" type="mesh" mesh="small_tet"/>
26+
</body>
27+
28+
<body name="tet3" pos="0.5 0 1.5">
29+
<joint name="root2" class="/" type="free"/>
30+
<geom name="tet3" contype="1" conaffinity="1" pos="0 0 0" type="mesh" mesh="small_tet"/>
31+
</body>
32+
33+
<body name="tet4" pos="0.1 0 2.0">
34+
<joint name="root3" class="/" type="free"/>
35+
<geom name="tet4" contype="1" conaffinity="1" pos="0 0 0" type="mesh" mesh="small_tet"/>
36+
</body>
37+
</worldbody>
38+
</mujoco>

genesis/engine/solvers/rigid/collider_decomp.py

Lines changed: 94 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -13,11 +13,20 @@
1313
from genesis.utils.misc import ti_field_to_torch
1414

1515
from .mpr_decomp import MPR
16+
from .gjk_decomp import GJK
17+
18+
from enum import IntEnum
1619

1720
if TYPE_CHECKING:
1821
from genesis.engine.solvers.rigid.rigid_solver_decomp import RigidSolver
1922

2023

24+
class CCD_ALGORITHM_CODE(IntEnum):
25+
MPR = 0
26+
MPR_SDF = 1
27+
GJK = 2
28+
29+
2130
@ti.func
2231
def rotaxis(vecin, i0, i1, i2, f0, f1, f2):
2332
vecres = ti.Vector([0.0, 0.0, 0.0], dt=gs.ti_float)
@@ -42,7 +51,18 @@ def __init__(self, rigid_solver: "RigidSolver"):
4251
self._solver = rigid_solver
4352
self._init_verts_connectivity()
4453
self._init_collision_fields()
54+
55+
# Identify the convex collision detection (ccd) algorithm
56+
if self._solver._options.use_gjk_collision:
57+
self.ccd_algorithm = CCD_ALGORITHM_CODE.GJK
58+
elif self._solver._enable_mujoco_compatibility:
59+
self.ccd_algorithm = CCD_ALGORITHM_CODE.MPR
60+
else:
61+
self.ccd_algorithm = CCD_ALGORITHM_CODE.MPR_SDF
62+
63+
# FIXME: MPR is necessary because it is used for terrain collision detection
4564
self._mpr = MPR(rigid_solver)
65+
self._gjk = GJK(rigid_solver) if self.ccd_algorithm == CCD_ALGORITHM_CODE.GJK else None
4666

4767
# multi contact perturbation and tolerance
4868
if self._solver._enable_mujoco_compatibility:
@@ -879,7 +899,7 @@ def _func_narrow_phase_convex_vs_convex(self):
879899
and self._solver.geoms_info[i_gb].type == gs.GEOM_TYPE.BOX
880900
)
881901
):
882-
self._func_mpr(i_ga, i_gb, i_b)
902+
self._func_convex_convex_contact(i_ga, i_gb, i_b)
883903

884904
@ti.kernel
885905
def _func_narrow_phase_convex_specializations(self):
@@ -897,7 +917,9 @@ def _func_narrow_phase_convex_specializations(self):
897917
and self._solver.geoms_info[i_gb].type == gs.GEOM_TYPE.BOX
898918
):
899919
if ti.static(sys.platform == "darwin"):
900-
self._func_mpr(i_ga, i_gb, i_b)
920+
# FIXME: It seems redundant, why don't we just call _func_plane_box_contact directly?
921+
# Anyway in this function, we will call _func_plane_box_contact.
922+
self._func_convex_convex_contact(i_ga, i_gb, i_b)
901923
else:
902924
self._func_plane_box_contact(i_ga, i_gb, i_b)
903925

@@ -1057,7 +1079,7 @@ def _func_plane_box_contact(self, i_ga, i_gb, i_b):
10571079
plane_dir = gu.ti_transform_by_quat(plane_dir, ga_state.quat)
10581080
normal = -plane_dir.normalized()
10591081

1060-
v1, _ = self._mpr.support_box(normal, i_gb, i_b)
1082+
v1, _ = self._mpr.support_field._func_support_box(normal, i_gb, i_b)
10611083
penetration = normal.dot(v1 - ga_state.pos)
10621084

10631085
if penetration > 0.0:
@@ -1173,7 +1195,7 @@ def _func_contact_orthogonals(self, i_ga, i_gb, normal, i_b):
11731195
return axis_0, axis_1
11741196

11751197
@ti.func
1176-
def _func_mpr(self, i_ga, i_gb, i_b):
1198+
def _func_convex_convex_contact(self, i_ga, i_gb, i_b):
11771199
if self._solver.geoms_info[i_ga].type > self._solver.geoms_info[i_gb].type:
11781200
i_ga, i_gb = i_gb, i_ga
11791201

@@ -1242,41 +1264,68 @@ def _func_mpr(self, i_ga, i_gb, i_b):
12421264
contact_pos = v1 - 0.5 * penetration * normal
12431265
is_col = penetration > 0
12441266
else:
1245-
# Try using MPR before anything else
1246-
is_mpr_updated = False
1247-
is_mpr_guess_direction_available = True
1248-
normal_ws = self.contact_cache[i_ga, i_gb, i_b].normal
1249-
for i_mpr in range(2):
1250-
if i_mpr == 1:
1251-
# Try without warm-start if no contact was detected using it.
1252-
# When penetration depth is very shallow, MPR may wrongly classify two geometries as not in
1253-
# contact while they actually are. This helps to improve contact persistence without increasing
1254-
# much the overall computational cost since the fallback should not be triggered very often.
1255-
is_mpr_guess_direction_available = (ti.abs(normal_ws) > gs.EPS).any()
1256-
if (i_detection == 0) and not is_col and is_mpr_guess_direction_available:
1257-
normal_ws = ti.Vector.zero(gs.ti_float, 3)
1258-
is_mpr_updated = False
1259-
1260-
if not is_mpr_updated:
1261-
is_col, normal, penetration, contact_pos = self._mpr.func_mpr_contact(
1262-
i_ga, i_gb, i_b, normal_ws
1263-
)
1264-
is_mpr_updated = True
1265-
1266-
# Fallback on SDF if collision is detected by MPR but no collision direction was cached but the
1267-
# initial penetration is already quite large, because the contact information provided by MPR
1268-
# may be unreliable in such a case.
1269-
# Here it is assumed that generic SDF is much slower than MPR, so it is faster in average
1270-
# to first make sure that the geometries are truly colliding and only after to run SDF if
1271-
# necessary. This would probably not be the case anymore if it was possible to rely on
1272-
# specialized SDF implementation for convex-convex collision detection in the first place.
1273-
if is_col and penetration > tolerance and not is_mpr_guess_direction_available:
1274-
# Note that SDF may detect different collision points depending on geometry ordering.
1275-
# Because of this, it is necessary to run it twice and take the contact information
1276-
# associated with the point of deepest penetration.
1277-
try_sdf = True
1278-
1279-
if ti.static(not self._solver._enable_mujoco_compatibility):
1267+
### MPR, MPR + SDF
1268+
if ti.static(self.ccd_algorithm != CCD_ALGORITHM_CODE.GJK):
1269+
# Try using MPR before anything else
1270+
is_mpr_updated = False
1271+
is_mpr_guess_direction_available = True
1272+
normal_ws = self.contact_cache[i_ga, i_gb, i_b].normal
1273+
for i_mpr in range(2):
1274+
if i_mpr == 1:
1275+
# Try without warm-start if no contact was detected using it.
1276+
# When penetration depth is very shallow, MPR may wrongly classify two geometries as not in
1277+
# contact while they actually are. This helps to improve contact persistence without increasing
1278+
# much the overall computational cost since the fallback should not be triggered very often.
1279+
is_mpr_guess_direction_available = (ti.abs(normal_ws) > gs.EPS).any()
1280+
if (i_detection == 0) and not is_col and is_mpr_guess_direction_available:
1281+
normal_ws = ti.Vector.zero(gs.ti_float, 3)
1282+
is_mpr_updated = False
1283+
1284+
if not is_mpr_updated:
1285+
is_col, normal, penetration, contact_pos = self._mpr.func_mpr_contact(
1286+
i_ga, i_gb, i_b, normal_ws
1287+
)
1288+
is_mpr_updated = True
1289+
1290+
# Fallback on SDF if collision is detected by MPR but no collision direction was cached but the
1291+
# initial penetration is already quite large, because the contact information provided by MPR
1292+
# may be unreliable in such a case.
1293+
# Here it is assumed that generic SDF is much slower than MPR, so it is faster in average
1294+
# to first make sure that the geometries are truly colliding and only after to run SDF if
1295+
# necessary. This would probably not be the case anymore if it was possible to rely on
1296+
# specialized SDF implementation for convex-convex collision detection in the first place.
1297+
if is_col and penetration > tolerance and not is_mpr_guess_direction_available:
1298+
# Note that SDF may detect different collision points depending on geometry ordering.
1299+
# Because of this, it is necessary to run it twice and take the contact information
1300+
# associated with the point of deepest penetration.
1301+
try_sdf = True
1302+
### GJK
1303+
elif ti.static(self.ccd_algorithm == CCD_ALGORITHM_CODE.GJK):
1304+
# If it was not the first detection, only detect single contact point.
1305+
self._gjk.func_gjk_contact(i_ga, i_gb, i_b, i_detection == 0)
1306+
1307+
is_col = self._gjk.is_col[i_b] == 1
1308+
penetration = self._gjk.penetration[i_b]
1309+
n_contacts = self._gjk.n_contacts[i_b]
1310+
1311+
if is_col:
1312+
if self._gjk.multi_contact_flag[i_b]:
1313+
# Used MuJoCo's multi-contact algorithm to find multiple contact points. Therefore,
1314+
# add the discovered contact points and stop multi-contact search.
1315+
for i_c in range(n_contacts):
1316+
if i_c >= self._n_contacts_per_pair:
1317+
# Ignore contact points if the number of contacts exceeds the limit.
1318+
break
1319+
contact_pos = self._gjk.contact_pos[i_b, i_c]
1320+
normal = self._gjk.normal[i_b, i_c]
1321+
self._func_add_contact(i_ga, i_gb, normal, contact_pos, penetration, i_b)
1322+
1323+
break
1324+
else:
1325+
contact_pos = self._gjk.contact_pos[i_b, 0]
1326+
normal = self._gjk.normal[i_b, 0]
1327+
1328+
if ti.static(self.ccd_algorithm == CCD_ALGORITHM_CODE.MPR_SDF):
12801329
if try_sdf:
12811330
is_col_a = False
12821331
is_col_b = False
@@ -1350,7 +1399,7 @@ def _func_mpr(self, i_ga, i_gb, i_b):
13501399
self.contact_cache[i_ga, i_gb, i_b].normal.fill(0.0)
13511400

13521401
elif multi_contact and is_col_0 > 0 and is_col > 0:
1353-
if ti.static(not self._solver._enable_mujoco_compatibility):
1402+
if ti.static(self.ccd_algorithm == CCD_ALGORITHM_CODE.MPR_SDF):
13541403
# 1. Project the contact point on both geometries
13551404
# 2. Revert the effect of small rotation
13561405
# 3. Update contact point
@@ -1387,6 +1436,11 @@ def _func_mpr(self, i_ga, i_gb, i_b):
13871436
# dynamics since zero-penetration contact points should not induce any force.
13881437
penetration = normal.dot(contact_point_b - contact_point_a)
13891438

1439+
elif ti.static(self.ccd_algorithm == CCD_ALGORITHM_CODE.GJK):
1440+
# Only change penetration to the initial one, because the normal vector could change abruptly
1441+
# under GJK-EPA as the nearest simplex is determined by discrete logic, unlike MPR.
1442+
penetration = penetration_0
1443+
13901444
# Discard contact point is repeated
13911445
repeated = False
13921446
for i_con in range(n_con):

0 commit comments

Comments
 (0)