Implementing a tetrahedra - ray intersection test

Question

I have a question regarding tetrahedra - ray intersections:

I tried to implement a tetrahedra-ray-intersection test which should return the index of the exit face. For the intersection, I followed this blog post, where scalar triple products were used for intersection testing:

http://realtimecollisiondetection.net/blog/?p=13

However, my code always returns the same face as exit face. After trying to find the solution myself without success, I would greatly appreciate any hints where the bug in my intersection test could lie.

Input parameters are the ray origin and direction, nodes, indices of the faces, indices of adjacent tetrahedra and index of last face. The indices of the exit face and exit tetrahedra are output values. ScTP computes the scalar triple product, and sameSign checks of the sign of three floats is the same. Here is the code, thanks for help:

__device__ void GetExitTet(float4 ray_o, float4 ray_d, float4* nodes, int32_t findex[4], int32_t adjtet[4], int32_t lface, int32_t &face, int32_t &tet)
{
face = 0;
tet = 0;
float4 q = ray_d;

    float4 v0 = make_float4(nodes[0].x, nodes[0].y, nodes[0].z, 0); // A
    float4 v1 = make_float4(nodes[1].x, nodes[1].y, nodes[1].z, 0); // B
    float4 v2 = make_float4(nodes[2].x, nodes[2].y, nodes[2].z, 0); // C
    float4 v3 = make_float4(nodes[3].x, nodes[3].y, nodes[3].z, 0); // D

    float4 p0 = v0 - ray_o;
    float4 p1 = v1 - ray_o;
    float4 p2 = v2 - ray_o;
    float4 p3 = v3 - ray_o;

    float u_3 = ScTP(q, p0, p1);
    float v_3 = ScTP(q, p1, p2);
    float w_3 = ScTP(q, p2, p0);

    float u_2 = ScTP(q, p1, p0);
    float v_2 = ScTP(q, p0, p3);
    float w_2 = ScTP(q, p3, p1);

    float u_1 = ScTP(q, p2, p3);
    float v_1 = ScTP(q, p3, p0);
    float w_1 = ScTP(q, p0, p2);

    float u_0 = ScTP(q, p3, p2);
    float v_0 = ScTP(q, p2, p1);
    float w_0 = ScTP(q, p1, p3);

    // ABC
    if (lface != findex[3]) { if (sameSign(u_3, v_3, w_3)) { face = findex[3]; tet = adjtet[3]; } }
    // BAD
    if (lface != findex[2]) { if (sameSign(u_2, v_2, w_2)) { face = findex[2]; tet = adjtet[2]; } }
    // CDA
    if (lface != findex[1]) { if (sameSign(u_1, v_1, w_1)) { face = findex[1]; tet = adjtet[1]; } }
    // DCB
    if (lface != findex[0]) { if (sameSign(u_0, v_0, w_0)) { face = findex[0]; tet = adjtet[0]; } }
    // No face hit
    // if (face == 0 && tet == 0) { printf("Error! No exit tet found. \n"); }
}

Answer 1

So, I came up with the solution myself, in case anyone needs it, here is it:

__device__ void GetExitTet(float4 ray_o, float4 ray_d, float4* nodes, int32_t findex[4], int32_t adjtet[4], int32_t lface, int32_t &face, int32_t &tet)
{
    face = 0;
    tet = 0;

    // http://realtimecollisiondetection.net/blog/?p=13

    // translate Ray to origin and vertices same as ray
    float4 q = ray_d;

    float4 v0 = make_float4(nodes[0].x, nodes[0].y, nodes[0].z, 0); // A
    float4 v1 = make_float4(nodes[1].x, nodes[1].y, nodes[1].z, 0); // B
    float4 v2 = make_float4(nodes[2].x, nodes[2].y, nodes[2].z, 0); // C
    float4 v3 = make_float4(nodes[3].x, nodes[3].y, nodes[3].z, 0); // D

    float4 p0 = v0 - ray_o;
    float4 p1 = v1 - ray_o;
    float4 p2 = v2 - ray_o;
    float4 p3 = v3 - ray_o;

    double QAB = ScTP(q, p0, p1); // A B
    double QBC = ScTP(q, p1, p2); // B C
    double QAC = ScTP(q, p0, p2); // A C
    double QAD = ScTP(q, p0, p3); // A D
    double QBD = ScTP(q, p1, p3); // B D
    double QCD = ScTP(q, p2, p3); // C D

    double sQAB = signf(QAB); // A B
    double sQBC = signf(QBC); // B C
    double sQAC = signf(QAC); // A C
    double sQAD = signf(QAD); // A D
    double sQBD = signf(QBD); // B D
    double sQCD = signf(QCD); // C D

    // ABC
    if (sQAB != 0 && sQAC !=0 && sQBC != 0) 
    { 
        if (sQAB < 0 && sQAC > 0 && sQBC < 0) { face = findex[3]; tet = adjtet[3]; } // exit face
    }
    // BAD
    if (sQAB != 0 && sQAD != 0 && sQBD != 0)
    {
        if (sQAB > 0 && sQAD < 0 && sQBD > 0) { face = findex[2]; tet = adjtet[2]; } // exit face
    }
    // CDA
    if (sQAD != 0 && sQAC != 0 && sQCD != 0)
    {
        if (sQAD > 0 && sQAC < 0 && sQCD < 0) { face = findex[1]; tet = adjtet[1]; } // exit face
    }
    // DCB
    if (sQBC != 0 && sQBD != 0 && sQCD != 0)
    {
        if (sQBC > 0 && sQBD < 0 && sQCD > 0) { face = findex[0]; tet = adjtet[0]; } // exit face
    }
    // No face hit
    // if (face == 0 && tet == 0) { printf("Error! No exit tet found. \n"); }
}