Line vs Triangle wrong computation

Question

I used the Möller-Trumbore algorithm. That seemed to work for a while but I found a bug:

const glm::vec3 avert0 = glm::vec3(-0.13f, 0.01f, 0.0f);
const glm::vec3 avert1 = glm::vec3(-0.01f, 0.01f, 0.0f);
const glm::vec3 avert2 = glm::vec3(-0.13f, 0.1f, 0.0f);

const glm::vec3 start = glm::vec3(-0.1f, 0.05f, -1.0f);
const glm::vec3 end = glm::vec3(-0.1f, 0.05f,    1.0f);

rayIntersectTriangle(start, end, avert0, avert1, avert2, position);

The above code does not collide ( u is < 0 ). However, if I change the Z value of the end vector from 1 to 30, it works:

 const glm::vec3 end = glm::vec3(-0.1f, 0.05f,    30.0f);

---

Here is my implementation:

bool Node::rayIntersectTriangle(glm::vec3 orig, glm::vec3 dest
              , glm::vec3 vert0, glm::vec3 vert1, glm::vec3 vert2
              , glm::vec3 &result)
{
const glm::vec3 edge1 = vert1 - vert0;
const glm::vec3 edge2 = vert2 - vert0;
const glm::vec3 pvec = glm::cross(dest, edge2);
const float det = glm::dot(edge1, pvec);

static const float Epsilon = std::numeric_limits<float>::epsilon();

/* 
   If det is below zero, we are back facing the triangle.
   If det is close to zero we are missing the triangle.
*/
if (det > -Epsilon && det < Epsilon)
  return false;

const float invDet = 1.0f / det;


const glm::vec3 tvec = orig - vert0;

result.x = glm::dot(tvec, pvec) * invDet;
if (result.x < 0.0f || result.x > 1.0f)
  return false;

const glm::vec3 qvec = glm::cross(tvec, edge1);

result.y = glm::dot(dest, qvec) * invDet;
if (result.y < 0.0f || result.x + result.y > 1.0f)
  return false;

result.z = glm::dot(edge2, qvec) * invDet;

return true;
}

Answer 1

I think the code only works for single sided tests. I have used a similar function before (adapted from the book: Real-Time Collision Detection by Christer Ericson) that does a double sided check:

template <typename T> int sgn(T val) {
    return (T(0) < val) - (val < T(0));
}

bool same_sign(const float a, const float b) {
    return sgn(a) == sgn(b);
}

float scalar_product(const glm::vec3& a, const glm::vec3& b, const glm::vec3& c) {
    return glm::dot(glm::cross(a, b), c);
}

bool IntersectLineTriangle(glm::vec3 p, glm::vec3 q, glm::vec3 a, glm::vec3 b, glm::vec3 c,
float &u, float &v, float &w)
{
    const glm::vec3 pq = q - p;
    const glm::vec3 pa = a - p;
    const glm::vec3 pb = b - p;
    const glm::vec3 pc = c - p;

    glm::vec3 m = glm::cross(pq, pc);
    u = glm::dot(pb, m); 
    v = -glm::dot(pa, m); 
    if (!same_sign(u, v)) return 0;
    w = scalar_product(pq, pb, pa);
    if (!same_sign(u, w)) return 0;

    const float denom = 1.0f / (u + v + w);
    u *= denom;
    v *= denom;
    w *= denom; // w = 1.0f - u - v;
    return true;
}

That returns a positive result for the inputs you give, although it may not be what you after.

EDIT1: possible explanation

I think the original function is giving a incorrect answer because the test works on the assumption that start->end passes through a triangle where the verts are defined counterclockwise. But the points above have the line segment going through the triangle from the other direction:

which effectively reverses the winding of the verts, from the point of view of the line segment. The double sided test takes this possibility into account and returns the proper result.

EDIT2: Original function should take direction, not end vector

After looking at the article I think I've found the error, I believe your original function should be taking the direction vector, rather than the end vector. It should look something like this:

bool rayIntersectTriangle(
    const glm::vec3& orig, const glm::vec3& dir, 
    const glm::vec3& vert0, const glm::vec3& vert1, const glm::vec3& vert2, 
    glm::vec3 &result, float& t)
{
    const glm::vec3 edge1 = vert1 - vert0;
    const glm::vec3 edge2 = vert2 - vert0;
    const glm::vec3 pvec = glm::cross(dir, edge2);
    const float det = glm::dot(edge1, pvec);

    static const float Epsilon = std::numeric_limits<float>::epsilon();

    /* 
       If det is below zero, we are back facing the triangle.
       If det is close to zero we are missing the triangle.
    */
    if (det > -Epsilon && det < Epsilon)
      return false;

    const float invDet = 1.0f / det;

    const glm::vec3 tvec = orig - vert0;

    result.x = glm::dot(tvec, pvec) * invDet;
    if (result.x < 0.0f || result.x > 1.0f)
      return false;

    const glm::vec3 qvec = glm::cross(tvec, edge1);

    result.y = glm::dot(dir, qvec) * invDet;
    if (result.y < 0.0f || (result.x + result.y) > 1.0f)
      return false; 

    result.z = 1.0f - result.x - result.y; // barycentric w = 1 - u - v
    t = glm::dot(edge2, qvec) * invDet; // ray intersection

    return true;
}

then call it like so:

rayIntersectTriangle(start, glm::normalize(end-start), avert0, avert1, avert2, position, t);