AlexanderFabisch / distance3d

Distance computation and collision detection in 3D.
https://alexanderfabisch.github.io/distance3d/
Other
51 stars 7 forks source link

Getting non reproducible collisions and penetration-depths #78

Open manuel-koch opened 1 year ago

manuel-koch commented 1 year ago

I have a huge set of meshes ( >5000 ) in my collision detection scene. If I execute gjk.gjk on those pairs of meshes with overlapping AABBs, I get ~227000 collisions. Filtering out the collisions that are of interest for my use-case, I still find false-positives with a penetration-depth of non-zero using epa.epa. Those false-positive collisions (sometimes) seem to be related to two meshes that are barely touching, see example scene.

two_cubes_barely_touching

If I re-run the collision again in a test script ( using the same vertex data of the two supposedly colliding meshes ) I get a collision with different penetration-depth. The difference is sometimes "huge", e.g. first time I get "4.0" after re-testing it I get "0.0". If they are indeed barely touching, I would expect a tiny penetration-depth in both cases.

Or if I shuffle the order of collision testing, I get different penetration-depth for the same pair of meshes. I.e. the following pseudo code of my own usage of distance3d yields different penetration-depth results for same mesh pairs that are in-collision:

  1. broad-phase, find potentially collision of meshes by overlapping AABBs of all meshes
    • e.g. this yields pair mesh indices [ (1,2), (5,8), (23,44), (5,15) ]
  2. find collisions using gjk.gjk and epa.epa between combination of AABB-overlapping-mesh-pair
    • i.e. I test mesh 1 vs 2, 5 vs 8, 23 vs 44, 5 vs 15
  3. shuffle list from 1.
    • e.g. this yields pair mesh indices [ (5,8), (1,2), (5,15), (23,44) ]
  4. redo 2. with the shuffled list of 3.
    • i.e. I test mesh 5 vs 8, 1 vs 2, 5 vs 15, 23 vs 44

results from 2. and 4. sometimes differ for the returned in-collision-flag and penetration-depth for the same pair of meshes !

Just swapping the two colliders or repeating the collisions yields different results for penetration-depth too.

def narrow_collision(collider_a, collider_b):
    distance, closest_point_a, closest_point_b, simplex = gjk.gjk(
        collider_a, collider_b
    )
    colliding = numpy.allclose(distance, 0)
    penetration_depth = 0.0
    if colliding:
        minimum_translation_vector, minkowski_faces, success = epa.epa(
            simplex,
            collider_a,
            collider_b,
            max_faces=512,
        )
        if success and not numpy.allclose(minimum_translation_vector, 0):
            penetration_depth = numpy.linalg.norm(minimum_translation_vector)
    return colliding, penetration_depth

colliding_one_two, penetration_depth_one_two = narrow_collision(
    my_collider_one, my_collider_two
)
colliding_two_one, penetration_depth_two_one = narrow_collision(
    my_collider_two, my_collider_one
)
# penetration_depth_one_two != penetration_depth_two_one

Summary: I see non reproducible results when

  1. using a different order of mesh-vs-mesh collisions
  2. or just swapping the meshes for a single mesh-vs-mesh collision
  3. or just repeating a single mesh-vs-mesh collision

e.g. examples of penetration-depths for repeated mesh-vs-mesh collisions I have seen ( the three cases are distinct pairs of meshes ):

Due to the vast amount of collisions I have to evaluate in my test case, it's hard to pin those collision result differences to just meshes that are barely touching. I'll have to further investigate on a mesh-by-mesh case.

Using:

AlexanderFabisch commented 1 year ago

Hi, thanks for reporting.

This might take some time to debug.

For now I'd propose to use mpr instead to compute the penetration depth.

manuel-koch commented 1 year ago

ok, I'll give it a try - thank you for the suggestion !

manuel-koch commented 1 year ago

@AlexanderFabisch Do I need to call gjk.gjk first to only call mpr in case of a collision ? Or can I just use mpr to find out both ( in-collision and penetration-depth ) in the same call ?

manuel-koch commented 1 year ago

Hm, even mpr.mpr_penetration returns different results, depending on the order of meshes. E.g. swapping the colliders in the mpr.mpr_penetration call and repeating it yields diffs like the following for distinct pairs of meshes ( not necessarily the same meshes as in the comments above ! ):

At least the category ( penetrating vs not-penetrating ) seems to be the same, even when swapping the colliders in the call.

AlexanderFabisch commented 1 year ago

Could you please create a minimal working example to reproduce this issue?

manuel-koch commented 1 year ago

Example script, resulting in a minor diff when swapping colliders:

A vs B: colliding=True penetration_depth=13.191347585870057
B vs A: colliding=True penetration_depth=13.013540835069607

generated by

from pyrr import Vector3
from trimesh import Trimesh, Scene

vertices_a = numpy.array(
    [Vector3([  40.     ,  157.98601, 1360.     ], dtype=numpy.float32), Vector3([  40.     ,  197.98601, 1360.     ], dtype=numpy.float32), Vector3([-2.6290081e-13,  1.9798601e+02,  1.3600000e+03], dtype=numpy.float32), Vector3([-2.6290081e-13,  1.5798601e+02,  1.3600000e+03], dtype=numpy.float32), Vector3([ 40.     , 157.98601,   0.     ], dtype=numpy.float32), Vector3([ 40.     , 197.98601,   0.     ], dtype=numpy.float32), Vector3([2.6290081e-13, 1.9798601e+02, 0.0000000e+00], dtype=numpy.float32), Vector3([2.6290081e-13, 1.5798601e+02, 0.0000000e+00], dtype=numpy.float32)],
    dtype=numpy.float64,
)
faces_a = numpy.array(
    [[0, 1, 3], [3, 1, 2], [5, 1, 0], [0, 4, 5], [0, 3, 7], [7, 4, 0], [2, 1, 6], [1, 5, 6], [6, 3, 2], [6, 7, 3], [6, 5, 4], [4, 7, 6]],
)
vertices_b = numpy.array(
    [Vector3([  24.        ,  154.98600769, 1340.6796875 ]), Vector3([  24.        ,  154.98600769, 1339.3203125 ]), Vector3([  24.        ,  170.98600769, 1339.3203125 ]), Vector3([  24.        ,  170.98600769, 1340.6796875 ]), Vector3([  19.20433235,  154.98600769, 1344.        ]), Vector3([  20.67968369,  154.98600769, 1344.        ]), Vector3([  20.67968369,  170.98600769, 1344.        ]), Vector3([  19.20433235,  170.98600769, 1344.        ]), Vector3([  16.        ,  154.98600769, 1340.7956543 ]), Vector3([  16.        ,  154.98600769, 1339.3203125 ]), Vector3([  16.54430389,  154.98600769, 1337.79882812]), Vector3([  17.7988472 ,  154.98600769, 1336.54431152]), Vector3([  19.32030869,  154.98600769, 1336.        ]), Vector3([  20.79566765,  154.98600769, 1336.        ]), Vector3([  22.2658329 ,  154.98600769, 1336.60900879]), Vector3([  23.45569611,  154.98600769, 1337.79882812]), Vector3([  23.45569992,  154.98600769, 1342.20117188]), Vector3([  22.20115662,  154.98600769, 1343.45568848]), Vector3([  17.7341671 ,  154.98600769, 1343.39099121]), Vector3([  16.60897255,  154.98600769, 1342.26586914]), Vector3([  19.32030869,  170.98600769, 1336.        ]), Vector3([  20.79566765,  170.98600769, 1336.        ]), Vector3([  16.        ,  170.98600769, 1339.3203125 ]), Vector3([  16.        ,  170.98600769, 1340.7956543 ]), Vector3([  16.60897255,  170.98600769, 1342.26586914]), Vector3([  17.7341671 ,  170.98600769, 1343.39099121]), Vector3([  22.20115662,  170.98600769, 1343.45568848]), Vector3([  23.45569992,  170.98600769, 1342.20117188]), Vector3([  23.45569611,  170.98600769, 1337.79882812]), Vector3([  22.2658329 ,  170.98600769, 1336.60900879]), Vector3([  17.7988472 ,  170.98600769, 1336.54431152]), Vector3([  16.54430389,  170.98600769, 1337.79882812]), Vector3([ 24.        , 154.98600769, -48.6796875 ]), Vector3([ 24.        , 170.98600769, -48.6796875 ]), Vector3([ 16.        , 154.98600769, -48.6796875 ]), Vector3([ 16.54430389, 154.98600769, -50.20117188]), Vector3([ 17.7988472 , 154.98600769, -51.45568848]), Vector3([ 19.32030869, 154.98600769, -52.        ]), Vector3([ 20.79566765, 154.98600769, -52.        ]), Vector3([ 22.2658329 , 154.98600769, -51.39099121]), Vector3([ 23.45569611, 154.98600769, -50.20117188]), Vector3([ 19.32030869, 170.98600769, -52.        ]), Vector3([ 20.79566765, 170.98600769, -52.        ]), Vector3([ 16.        , 170.98600769, -48.6796875 ]), Vector3([ 23.45569611, 170.98600769, -50.20117188]), Vector3([ 22.2658329 , 170.98600769, -51.39099121]), Vector3([ 17.7988472 , 170.98600769, -51.45568848]), Vector3([ 16.54430389, 170.98600769, -50.20117188])],
    dtype=numpy.float64,
)
faces_b = numpy.array(
    [[0, 1, 3], [3, 1, 2], [3, 16, 0], [3, 27, 16], [26, 5, 17], [6, 5, 26], [17, 16, 26], [16, 27, 26], [9, 8, 23], [22, 9, 23], [23, 8, 19], [19, 24, 23], [18, 8, 9], [12, 13, 18], [18, 11, 12], [18, 5, 4], [0, 16, 18], [14, 15, 18], [18, 16, 17], [18, 1, 0], [9, 10, 18], [17, 5, 18], [15, 1, 18], [19, 8, 18], [18, 13, 14], [10, 11, 18], [18, 24, 19], [18, 25, 24], [4, 5, 7], [7, 5, 6], [7, 18, 4], [25, 18, 7], [2, 28, 30], [6, 26, 30], [30, 31, 22], [30, 21, 20], [30, 3, 2], [22, 23, 30], [30, 29, 21], [30, 7, 6], [27, 3, 30], [30, 26, 27], [30, 23, 24], [24, 25, 30], [28, 29, 30], [25, 7, 30], [41, 42, 38], [37, 41, 38], [44, 32, 40], [33, 32, 44], [47, 35, 34], [47, 34, 43], [39, 38, 45], [45, 38, 42], [45, 40, 39], [45, 44, 40], [46, 41, 37], [37, 36, 46], [35, 47, 46], [46, 36, 35], [32, 33, 1], [1, 33, 2], [34, 35, 9], [9, 35, 10], [35, 36, 10], [10, 36, 11], [36, 37, 11], [11, 37, 12], [37, 38, 12], [12, 38, 13], [38, 39, 13], [13, 39, 14], [40, 32, 15], [15, 32, 1], [39, 40, 14], [14, 40, 15], [42, 41, 21], [21, 41, 20], [43, 34, 22], [22, 34, 9], [33, 44, 2], [2, 44, 28], [45, 42, 29], [29, 42, 21], [44, 45, 28], [28, 45, 29], [41, 46, 20], [20, 46, 30], [47, 43, 31], [31, 43, 22], [46, 47, 30], [30, 47, 31]],
)

collider_a = colliders.ConvexHullVertices(vertices_a)
collider_b = colliders.ConvexHullVertices(vertices_b)
colliding, penetration_depth, penetration_dir, contact_position = mpr.mpr_penetration(
    collider_a, collider_b
)
print(f"A vs B: colliding={colliding} penetration_depth={penetration_depth}")

colliding, penetration_depth, penetration_dir, contact_position = mpr.mpr_penetration(
    collider_b, collider_a
)
print(f"B vs A: colliding={colliding} penetration_depth={penetration_depth}")

scene = Scene()
mesh_a = Trimesh(vertices=vertices_a, faces=faces_a)
mesh_a.visual.face_colors = (255, 0, 0, 80)
scene.add_geometry(mesh_a)
mesh_b = Trimesh(vertices=vertices_b, faces=faces_b)
mesh_b.visual.face_colors = (0, 255, 0, 80)
scene.add_geometry(mesh_b)
scene.show()

The scene looks like

example_1_solid example_1_wireframe

===============================================

Example script, resulting in a major diff when swapping colliders:

A vs B: colliding=True penetration_depth=4.50000000000012
B vs A: colliding=True penetration_depth=14.411558781415797

generated by

    from pyrr import Vector3
    from trimesh import Trimesh, Scene

    vertices_a = numpy.array(
        [Vector3([  40.     ,  157.98601, 1360.     ], dtype=numpy.float32), Vector3([  40.     ,  197.98601, 1360.     ], dtype=numpy.float32), Vector3([-2.6290081e-13,  1.9798601e+02,  1.3600000e+03], dtype=numpy.float32), Vector3([-2.6290081e-13,  1.5798601e+02,  1.3600000e+03], dtype=numpy.float32), Vector3([ 40.     , 157.98601,   0.     ], dtype=numpy.float32), Vector3([ 40.     , 197.98601,   0.     ], dtype=numpy.float32), Vector3([2.6290081e-13, 1.9798601e+02, 0.0000000e+00], dtype=numpy.float32), Vector3([2.6290081e-13, 1.5798601e+02, 0.0000000e+00], dtype=numpy.float32)],
        dtype=numpy.float64,
    )
    faces_a = numpy.array(
        [[0, 1, 3], [3, 1, 2], [5, 1, 0], [0, 4, 5], [0, 3, 7], [7, 4, 0], [2, 1, 6], [1, 5, 6], [6, 3, 2], [6, 7, 3], [6, 5, 4], [4, 7, 6]],
    )
    vertices_b = numpy.array(
        [Vector3([  4.5       , 328.48599243, 989.49023438]), Vector3([  4.5       , 328.48599243, 990.50976562]), Vector3([-29.5       , 328.48599243, 990.50976562]), Vector3([-29.5       , 328.48599243, 989.49023438]), Vector3([  4.5       , 332.08276367, 987.        ]), Vector3([  4.5       , 330.97622681, 987.        ]), Vector3([-29.5       , 330.97622681, 987.        ]), Vector3([-29.5       , 332.08276367, 987.        ]), Vector3([  4.5       , 334.48599243, 989.40325928]), Vector3([  4.5       , 334.48599243, 990.50976562]), Vector3([  4.5       , 334.07775879, 991.65087891]), Vector3([  4.5       , 333.13687134, 992.59179688]), Vector3([  4.5       , 331.99575806, 993.        ]), Vector3([  4.5       , 330.88925171, 993.        ]), Vector3([  4.5       , 329.78662109, 992.54327393]), Vector3([  4.5       , 328.89422607, 991.65087891]), Vector3([  4.5       , 328.89422607, 988.34912109]), Vector3([  4.5       , 329.83514404, 987.40820312]), Vector3([  4.5       , 333.18536377, 987.45672607]), Vector3([  4.5       , 334.02926636, 988.30059814]), Vector3([-29.5       , 331.99575806, 993.        ]), Vector3([-29.5       , 330.88925171, 993.        ]), Vector3([-29.5       , 334.48599243, 990.50976562]), Vector3([-29.5       , 334.48599243, 989.40325928]), Vector3([-29.5       , 334.02926636, 988.30059814]), Vector3([-29.5       , 333.18536377, 987.45672607]), Vector3([-29.5       , 329.83514404, 987.40820312]), Vector3([-29.5       , 328.89422607, 988.34912109]), Vector3([-29.5       , 328.89422607, 991.65087891]), Vector3([-29.5       , 329.78662109, 992.54327393]), Vector3([-29.5       , 333.13687134, 992.59179688]), Vector3([-29.5       , 334.07775879, 991.65087891]), Vector3([  4.5       , 182.48599243, 989.49023438]), Vector3([  4.5       , 182.48599243, 990.50976562]), Vector3([-29.5       , 182.48599243, 990.50976562]), Vector3([-29.5       , 182.48599243, 989.49023438]), Vector3([  4.5       , 184.97622681, 987.        ]), Vector3([-29.5       , 184.97622681, 987.        ]), Vector3([  4.5       , 184.88925171, 993.        ]), Vector3([  4.5       , 183.78662109, 992.54327393]), Vector3([  4.5       , 182.89422607, 991.65087891]), Vector3([  4.5       , 182.89422607, 988.34912109]), Vector3([  4.5       , 183.83514404, 987.40820312]), Vector3([-29.5       , 184.88925171, 993.        ]), Vector3([-29.5       , 183.83514404, 987.40820312]), Vector3([-29.5       , 182.89422607, 988.34912109]), Vector3([-29.5       , 182.89422607, 991.65087891]), Vector3([-29.5       , 183.78662109, 992.54327393])],
        dtype=numpy.float64,
    )
    faces_b = numpy.array(
        [[20, 21, 13], [12, 20, 13], [31, 10, 9], [31, 9, 22], [9, 8, 23], [22, 9, 23], [23, 8, 19], [19, 24, 23], [18, 8, 9], [12, 13, 18], [18, 11, 12], [18, 5, 4], [0, 16, 18], [14, 15, 18], [18, 16, 17], [18, 1, 0], [9, 10, 18], [17, 5, 18], [15, 1, 18], [19, 8, 18], [18, 13, 14], [10, 11, 18], [18, 24, 19], [18, 25, 24], [4, 5, 7], [7, 5, 6], [7, 18, 4], [25, 18, 7], [2, 28, 30], [6, 26, 30], [30, 31, 22], [30, 21, 20], [30, 3, 2], [22, 23, 30], [30, 29, 21], [30, 7, 6], [27, 3, 30], [30, 26, 27], [30, 23, 24], [24, 25, 30], [28, 29, 30], [25, 7, 30], [30, 20, 12], [12, 11, 30], [10, 31, 30], [30, 11, 10], [46, 33, 40], [34, 33, 46], [32, 33, 35], [35, 33, 34], [35, 41, 32], [35, 45, 41], [44, 36, 42], [37, 36, 44], [42, 41, 44], [41, 45, 44], [39, 38, 47], [47, 38, 43], [47, 40, 39], [47, 46, 40], [33, 32, 1], [1, 32, 0], [35, 34, 3], [3, 34, 2], [36, 37, 5], [5, 37, 6], [38, 39, 13], [13, 39, 14], [40, 33, 15], [15, 33, 1], [39, 40, 14], [14, 40, 15], [32, 41, 0], [0, 41, 16], [42, 36, 17], [17, 36, 5], [41, 42, 16], [16, 42, 17], [43, 38, 21], [21, 38, 13], [37, 44, 6], [6, 44, 26], [45, 35, 27], [27, 35, 3], [44, 45, 26], [26, 45, 27], [34, 46, 2], [2, 46, 28], [47, 43, 29], [29, 43, 21], [46, 47, 28], [28, 47, 29]],
    )

    collider_a = colliders.ConvexHullVertices(vertices_a)
    collider_b = colliders.ConvexHullVertices(vertices_b)
    colliding, penetration_depth, penetration_dir, contact_position = mpr.mpr_penetration(
        collider_a, collider_b
    )
    print(f"A vs B: colliding={colliding} penetration_depth={penetration_depth}")

    colliding, penetration_depth, penetration_dir, contact_position = mpr.mpr_penetration(
        collider_b, collider_a
    )
    print(f"B vs A: colliding={colliding} penetration_depth={penetration_depth}")

    scene = Scene()
    mesh_a = Trimesh(vertices=vertices_a, faces=faces_a)
    mesh_a.visual.face_colors = (255, 0, 0, 80)
    scene.add_geometry(mesh_a)
    mesh_b = Trimesh(vertices=vertices_b, faces=faces_b)
    mesh_b.visual.face_colors = (0, 255, 0, 80)
    scene.add_geometry(mesh_b)
    scene.show()

The scene looks like

example_2_solid example_2_wireframe
manuel-koch commented 1 year ago

Example script, resulting in a major diff when swapping colliders:

A vs B: colliding=True penetration_depth=11.800003050429694
B vs A: colliding=True penetration_depth=198.80704230984915

generated by

    from pyrr import Vector3
    from trimesh import Trimesh, Scene

    vertices_a = numpy.array(
        [Vector3([-670.     ,  157.98601,   40.     ], dtype=numpy.float32), Vector3([-670.     ,  197.98601,   40.     ], dtype=numpy.float32), Vector3([-6.7000000e+02,  1.9798601e+02,  2.8421709e-14], dtype=numpy.float32), Vector3([-6.7000000e+02,  1.5798601e+02,  1.4210855e-14], dtype=numpy.float32), Vector3([  0.     , 157.98601,  40.     ], dtype=numpy.float32), Vector3([  0.     , 197.98601,  40.     ], dtype=numpy.float32), Vector3([ 0.0000000e+00,  1.9798601e+02, -1.4210855e-14], dtype=numpy.float32), Vector3([ 0.0000000e+00,  1.5798601e+02, -2.8421709e-14], dtype=numpy.float32)],
        dtype=numpy.float64,
    )
    faces_a = numpy.array(
        [[0, 1, 3], [3, 1, 2], [5, 1, 0], [0, 4, 5], [0, 3, 7], [7, 4, 0], [2, 1, 6], [1, 5, 6], [6, 3, 2], [6, 7, 3], [6, 5, 4], [4, 7, 6]],
    )
    vertices_b = numpy.array(
        [Vector3([-197.        ,  193.48600769,   24.        ]), Vector3([-197.        ,  186.18600464,   24.        ]), Vector3([-197.        ,  186.18600464,   16.        ]), Vector3([-197.        ,  193.48600769,   16.        ]), Vector3([-210.        ,  193.48600769,   24.        ]), Vector3([-210.        ,  186.18600464,   24.        ]), Vector3([-210.        ,  186.18600464,   16.        ]), Vector3([-210.        ,  193.48600769,   16.        ]), Vector3([ 19.        , 193.48600769,  24.        ]), Vector3([ 19.        , 186.18600464,  24.        ]), Vector3([ 19.        , 186.18600464,  16.        ]), Vector3([ 19.        , 193.48600769,  16.        ])],
        dtype=numpy.float64,
    )
    faces_b = numpy.array(
        [[5, 1, 0], [0, 4, 5], [0, 3, 7], [7, 4, 0], [2, 1, 6], [1, 5, 6], [6, 3, 2], [6, 7, 3], [6, 5, 4], [4, 7, 6], [8, 9, 11], [11, 9, 10], [9, 8, 1], [1, 8, 0], [10, 9, 2], [2, 9, 1], [8, 11, 0], [0, 11, 3], [11, 10, 3], [3, 10, 2]],
    )

    collider_a = colliders.ConvexHullVertices(vertices_a)
    collider_b = colliders.ConvexHullVertices(vertices_b)
    colliding, penetration_depth, penetration_dir, contact_position = mpr.mpr_penetration(
        collider_a, collider_b
    )
    print(f"A vs B: colliding={colliding} penetration_depth={penetration_depth}")

    colliding, penetration_depth, penetration_dir, contact_position = mpr.mpr_penetration(
        collider_b, collider_a
    )
    print(f"B vs A: colliding={colliding} penetration_depth={penetration_depth}")

    scene = Scene()
    mesh_a = Trimesh(vertices=vertices_a, faces=faces_a)
    mesh_a.visual.face_colors = (255, 0, 0, 80)
    scene.add_geometry(mesh_a)
    mesh_b = Trimesh(vertices=vertices_b, faces=faces_b)
    mesh_b.visual.face_colors = (0, 255, 0, 80)
    scene.add_geometry(mesh_b)
    scene.show()

The scene looks like

example_3_solid example_3_wireframe
AlexanderFabisch commented 1 year ago

I looked into this case:

Example script, resulting in a major diff when swapping colliders:

A vs B: colliding=True penetration_depth=11.800003050429694
B vs A: colliding=True penetration_depth=198.80704230984915

The problem is that both solutions are actually correct. I am not aware of any guarantee for the penetration depth along the found direction to be minimal or unique. The amount of publications on MPR is rather limited. I guess the situation with EPA is the same although I didn't check.

If you are interested, I would suggest to look a bit more into the literature to find out more about penetration depth computation. I'm also curious if there are better algorithms or if I am missing anything here.

Another possible solution using distance 3d to find a unique penetration direction would be the hydroelastic contact model. However, it might give you a unique direction, but not a depth. It's also not trivial to compute the required tetrahedral mesh for arbitrary triangular meshes. If you know the underlying parameters of the primitive shape (box, cylinder, etc.) it would be possible to compute it though.

manuel-koch commented 1 year ago

Thanks for looking into the issue and providing some background info.

manuel-koch commented 1 year ago

For my use-case I will workaround the issue by doing collision-test twice ( a vs b and b vs a ) and use the minimal penetration-depth ( if available ) of the two results.

AlexanderFabisch commented 1 year ago

That's a good solution I guess. :)