Collision algorithms are not numerically robust
After many tweaks of tolerances in new geometric library, one can everytime find even worse conditioned simplices such that the algorithms stop working properly. Indeed, in the example below meshes are of very poor quality. So the question is: should these issues be handled more robustly (I can think of exact predicates in CollisionDetection::collides_(triangle|tetrahedron)_point
) or has a user to reconcile to strange failures on poor meshes?
from dolfin import *
def test_collision_robustness(aspect, y, step=0.00001):
_int_limit = 2**32 - 1
nx = 10
ny = int(aspect*nx)
mesh = UnitSquareMesh(nx, ny, 'crossed')
bb = mesh.bounding_box_tree()
print "Testing collision of points (x = %e:%e:%e, y = %e)" \
" with UnitSquareMesh %ix%i crossed:" \
% (0.0, 1.0, step, y, nx, ny)
print " min radius ratio", MeshQuality().radius_ratio_min_max(mesh)[0]
x = 0.0
while x <= 1.0:
p = Point(x, y)
c = bb.compute_first_entity_collision(p)
if c >= _int_limit:
print ' No collision at point', p.str()
closest = Vertex(mesh, 0).point()
distance = p.distance(closest)
for v in vertices(mesh):
if p.distance(v.point()) <= distance:
closest = v.point()
distance = p.distance(closest)
print ' Closest mesh vertex at point', closest.str(), ', distance', distance
return p, closest
x += step
print " OK"
if __name__=='__main__':
test_collision_robustness( 100, 1e-14)
test_collision_robustness( 40, 1e-03)
test_collision_robustness( 10, 1e-16, 1e-7)
test_collision_robustness(4.43, 1e-17, 4.03e-6)
test_collision_robustness( 100, 0.5 + 1e-14)
test_collision_robustness( 40, 0.5 + 1e-17, 1e-6)
test_collision_robustness( 10, 0.5 + 1e-16, 1e-7)
test_collision_robustness(4.43, 0.5 + 1e-17, 4.03e-6)
All the four tests fail on my machine.
Comments (10)
-
-
reporter - edited description
Added test cases failing in the interior of domain.
-
reporter @logg Ok, I agree with you. You're on thin ice when testing a point near facet. But we should guarantee that point is 'at least on one side of facet' by construction of an algorithm. Try this test
test_collision_robustness( 100, 0.5 + 1e-14) test_collision_robustness( 40, 0.5 + 1e-17, 1e-6) test_collision_robustness( 10, 0.5 + 1e-16, 1e-7) test_collision_robustness(4.43, 0.5 + 1e-17, 4.03e-6)
It shows that points which are 'very much' inside the domain but do not belong to any cell. We should guarantee that they belong to at least one. I will think about it.
-
I think it should be fairly easy to fix the algorithm to use the first method described in http://www.blackpawn.com/texts/pointinpoly/ which shouldn't allow a point to slip between two facets. e.g. something like this... though it probably needs some thought about correctness on manifolds.
bool CollisionDetection::collides_triangle_point(const Point& p0, const Point& p1, const Point& p2, const Point &point) { const Point norm = (p1 - p0).cross(p0 - p2); Point p = point - p0; double t1 = norm.dot(p.cross(p0 - p2)); if (t1 < 0) return false; p = point - p1; double t2 = norm.dot(p.cross(p1 - p0)); if (t2 < 0) return false; p = point - p2; double t3 = norm.dot(p.cross(p2 - p1)); if (t3 < 0) return false; return true; }
-
See proposed fix in pull request #186.
-
- changed milestone to 1.5
-
assigned issue to
-
- changed status to resolved
Fixed in branch chris/in-triangle, now merged.
-
There seems to still be issues related to this, the following fails for me:
from dolfin import * mesh = Mesh() editor = MeshEditor() editor.open(mesh, 1, 3) editor.init_vertices(2) editor.init_cells(1) editor.add_vertex(0, 41.06309891, 63.74219894, 68.10320282) editor.add_vertex(1, 41.45830154, 62.61560059, 66.43019867) editor.add_cell(0,0,1) editor.close() cell = Cell(mesh, 0) assert cell.contains(cell.midpoint())
-
I haven't looked at 1D, and manifold meshes are always more difficult to deal with. It may be a case of adjusting some tolerances in this example.
-
- removed milestone
Removing milestone: 1.5 (automated comment)
- Log in to comment
These tests all boil down to testing whether a point located on the edge of a triangle is inside the triangle. The test in question is implemented in the function
collides_triangle_point
indolfin/geometry/CollisionDetection.cpp
:This is a standard test (see for example "Real-Time Collision Detection" by Christer Ericson) and it is bound to fail in "extreme" cases as a result of round-off error.
Perhaps it can be tweaked to be more robust, but I don't see how it can be made completely robust without resorting to exact numerics. The applications I have in mind for collision detection are those where if one should miss the intersections of two triangles that collide only at a single point (or have a vanishingly small area of the overlap), then it doesn't matter for the numerical method.
Suggestions for an alternative more robust implementation of collides_triangle_point are welcome. Perhaps we could consider adding a more expensive version which would be optional (based on some parameter), but it's not of high priority for me to delve into.