Ray-Triangle Intersection C++

Question

I am testing if a ray intersects a triangle. For the time being, I'm using the following code to test if there is an intersection between a triangle, and the ray from origin to the midpoint of the triangle:

Ray<float> *ray = new Ray<float>(Vec3<float>(0), chosenTriangle->GetTriangleMidpoint()); 

Along side is the Vec3 object which I'm using to store the vector operations:

template<typename T>
class Vec3
{
public:
    T x, y, z;
    Vec3() : x(T(0)), y(T(0)), z(T(0)) { }
    Vec3(T xx) : x(xx), y(xx), z(xx) { }

    Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {}
    Vec3& normalize()
    {
        T nor2 = length2();
        if (nor2 > 0) {
            T invNor = 1 / sqrt(nor2);
            x *= invNor, y *= invNor, z *= invNor;
        }
        return *this;
    }

    Vec3<T> operator * (const T &f) const { return Vec3<T>(x * f, y * f, z * f); }
    Vec3<T> operator * (const Vec3<T> &v) const { return Vec3<T>(x * v.x, y * v.y, z * v.z); }
    T dot(const Vec3<T> &v) const { return x * v.x + y * v.y + z * v.z; }
    Vec3<T> operator - (const Vec3<T> &v) const { return Vec3<T>(x - v.x, y - v.y, z - v.z); }
    Vec3<T> operator + (const Vec3<T> &v) const { return Vec3<T>(x + v.x, y + v.y, z + v.z); }
    bool operator == (const Vec3<T> &v) { return x == v.x && y == v.y && z == v.z; }
    Vec3<T> operator - () const { return Vec3<T>(-x, -y, -z); }
    T length2() const { return x * x + y * y + z * z; }
    T length() const { return sqrt(length2()); }
    Vec3<T> CrossProduct(Vec3<T> other)
    { 
        return Vec3<T>(y*other.z - other.y*z, x*other.z - z*other.x, x*other.y - y*other.x); 
    }
    friend std::ostream & operator << (std::ostream &os, const Vec3<T> &v)
    {
        os << "[" << v.x << " " << v.y << " " << v.z << "]";
        return os;
    }

The chosen triangle and the ray have the following values, where vertA, vertB and vertC are the vertices of the triangle and are found in an object which represents a triangle.

The code which computes if there is an intersection between a specified ray and an intersection is the following. This code is found inside the triangle object method where vertA, vertB and vertC are global variables.

bool CheckRayIntersection(Vec3<T> &o, Vec3<T> &d)
{
    Vec3<T> e1 = vertB - vertA;
    Vec3<T> e2 = vertC - vertA;
    Vec3<T> p = d.CrossProduct(e2);
    T a = e1.dot(p);

    if(a == 0)
        return false;

    float f = 1.0f/a;

    Vec3<T> s = o - vertA;
    T u = f * s.dot(p);
    if(u < 0.0f || u > 1.0f)
        return false;

    Vec3<T> q = s.CrossProduct(e1);
    T v = f * d.dot(q);

    if(v < 0.0f || u+v > 1.0f)
        return false;

    T t = f * e2.dot(q);

    return (t >= 0);
    
}

I still get a false returned from the function, but I'm presuming it should return a true since a vector passing through the midpoint of the triangle should intersect the triangle at the midpoint. Can anybody enlighten me what's wrong in my code? Or is the returned false actually correct?

Answer 1

There is the (experimental) glm::intersectRayTriangle() function in GLM's GLM_GTX_intersect extension, you'll need to #define GLM_ENABLE_EXPERIMENTAL and #include <glm/gtx/intersect.hpp> to access it.

Answer 2

Direct translation of MichaelCMS's answer for use with glm.

// must normalize direction of ray
bool intersectRayTri(Tri& tri, glm::vec3 o, glm::vec3 n) {
    glm::vec3 e1, e2, pvec, qvec, tvec;

    e1 = tri.v2 - tri.v1;
    e2 = tri.v3 - tri.v1;
    pvec = glm::cross(n, e2);
    
    n = glm::normalize(n);
    //NORMALIZE(pvec);
    float det = glm::dot(pvec, e1);

    if (det != 0)
    {
        float invDet = 1.0f / det;
        tvec = o - tri.v1;
        // NORMALIZE(tvec);
        float u = invDet * glm::dot(tvec, pvec);
        if (u < 0.0f || u > 1.0f)
        {

            return false;
        }
        qvec = glm::cross(tvec, e1);
        // NORMALIZE(qvec);
        float v = invDet * glm::dot(qvec, n);
        if (v < 0.0f || u + v > 1.0f)
        {

            return false;
        }
    }
    else return false;
    return true; // det != 0 and all tests for false intersection fail
}

Answer 3

With your data, I managed to get consistent results by having the ray direction normalized (this is the only apparent change in the code).

Here is the code implementation (I used the paper as reference, and it's not very optimized) :

struct quickVect
{

   float x,y,z;
   float l;
 };

#define DOT(v1,v2) (v1.x*v2.x + v1.y*v2.y+v1.z*v2.z)
#define CROSS(rez,v1,v2) \
rez.x  = v1.y*v2.z - v1.z*v2.y; \
rez.y  = v1.z*v2.x - v1.x*v2.z; \
rez.z  = v1.x*v2.y - v1.y*v2.x;

#define SUB(rez,v1,v2) \
rez.x = v1.x-v2.x; \
rez.y = v1.y-v2.y; \
rez.z = v1.z-v2.z;


#define LENGTH(v) (sqrtf(v.x* v.x + v.y*v.y + v.z*v.z))

#define NORMALIZE(v) \
v.l = LENGTH(v); \
v.x = v.x / v.l; \
v.y = v.y / v.l; \
v.z = v.z / v.l;

#define EPSILON 0.000001f

//#define TEST_CULL

bool testIntersection(quickVect& v1, quickVect& v2, quickVect& v3, quickVect& orig,quickVect& dir)
{
 quickVect e1,e2,pvec,qvec,tvec;

 SUB(e1,v2,v1);
 SUB(e2,v3,v1);

 CROSS(pvec,dir,e2);

 NORMALIZE(dir);
 //NORMALIZE(pvec);
 float det = DOT(pvec,e1);
#ifdef TEST_CULL
if (det <EPSILON)
{

    return false;
}
SUB(tvec,orig,v1);
float u = DOT(tvec,pvec);
if (u < 0.0 || u > det)
{

    return false;
}
CROSS(qvec,tvec,e1);
float v = DOT(dir,qvec);
if (v < 0.0f || v + u > det)
{

    return false;
}
#else
 if (det < EPSILON && det > -EPSILON )
 {

     return false;
 }

 float invDet = 1.0f / det;
 SUB(tvec,orig,v1);
// NORMALIZE(tvec);
 float u = invDet * DOT(tvec,pvec);
 if (u <0.0f || u > 1.0f)
 {

     return false;
 }
 CROSS(qvec,tvec,e1);
// NORMALIZE(qvec);
 float v = invDet* DOT(qvec,dir);
 if (v < 0.0f || u+v > 1.0f)
 {

     return false;
 }
#endif
 return true;
}