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?

Aldredge

Where did vertA, vertB and vertC came from? They are not within function parameters. Meaning those are global variables or you didn't give complete function. – Subtangent

@Subtangent didnt give the complete code, they are just the coordinates of the triangle , will update – Aldredge

It seems you are implementing the Möller–Trumbore algorithm. My gut tells me you got one of your Vec3's methods wrong, since your implementation has the shape of the one I referred to – Neill

You should normalize e1 and e2 before doing cross product. – Objection

The paper referenced by Rerito says that d should be normalized. – Dissident

@Rerito: GOod to know. Wanted to write him an "unoptimized" step-by-step collision detection, but gave up because it became too long. – Subtangent

If you can't make algorithm Rerito linked, you could try treating triangle as 4 planes. It'll be less efficieant and longer, but easier to understand. It goes like this: 1. Calculate ray-plane intersection point for plane defined by triangle. 2. Compare against triangle's sides. Each side is defined as plane with "normal" pointing "inside" of the triangle. To have collision you need to have dot((point - planeOrigin), planeNormal) > 0.0 for all 3 planes created by triangle sides. – Subtangent

At what condition does it exit ? – Objection

@AdrianDeBarro Considering your test case, could you tell us which test yields the return false;? – Neill

@Neill its the second if statment i.e. if(u < 0.0f || u > 1.0f) return false; – Aldredge

@Objection will add the normalizations to the code and will check the vec3 methods just in case :) – Aldredge

Given your test case, you should get u == v == 1f/3f and t == 1 – Neill

@AdrianDeBarro after reading the paper, it seems that it doesn't use normalized directions (which I personally find weird). Is your code exiting because u > 1.0f ? (this is the usual fall if you use dot on un normalized directions ). I see however that a determinat is computed for normalization reasons in the algorithm. – Objection

Here's the paper in case you want exact implementations on the vector methods : cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/… – Objection

@Objection its failing at u < 0 whgich is -122. something – Aldredge

@Neill added the vec3 methods just in case but not seeing anything that its not normal in it. – Aldredge

@AdrianDeBarro : really off topic, but check out this raytracing tutorial : flipcode.com/archives/… . Check out the ray / triangle interesection implemented there (it has code as well). You can at least check if your result is consistend with another implementation. – Objection

@Objection good idea thanks will work on it – Aldredge

Answer 1

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

Objection

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
}

Mosley