Collision algorithms are not numerically robust

[Jan Blechta]

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.

[Anders Logg (Chalmers)]

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 in dolfin/geometry/CollisionDetection.cpp:

// Compute vectors
const Point v1 = p1 - p0;
const Point v2 = p2 - p0;
const Point v = point - p0;

// Compute entries of linear system
const double a11 = v1.dot(v1);
const double a12 = v1.dot(v2);
const double a22 = v2.dot(v2);
const double b1 = v.dot(v1);
const double b2 = v.dot(v2);

// Solve linear system
const double inv_det = 1.0 / (a11*a22 - a12*a12);
const double x1 = inv_det*( a22*b1 - a12*b2);
const double x2 = inv_det*(-a12*b1 + a11*b2);

// Tolerance for numeric test (using vector v1)
const double dx = std::abs(v1.x());
const double dy = std::abs(v1.y());
const double eps = std::max(DOLFIN_EPS_LARGE,    DOLFIN_EPS_LARGE*std::max(dx, dy));

// Check if point is inside
return x1 >= -eps && x2 >= -eps && x1 + x2 <= 1.0 + eps;

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.

[Jan Blechta]

Added test cases failing in the interior of domain.

[Jan Blechta]

@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.

[Chris Richardson]

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;
}

[Chris Richardson]

See proposed fix in pull request #186.

[Chris Richardson]

Fixed in branch chris/in-triangle, now merged.

[Øyvind Evju]

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())

[Chris Richardson]

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.

[Prof Garth Wells]

Removing milestone: 1.5 (automated comment)