射线与三角面元求交点

#include <CSTDIO>
#include <cmath>

#define EPSILON 0.000001
/*#define CROSS(dest, v1, v2) \
    dest[0] = v1[1]*v2[2] - v1[2]*v2[1];\
    dest[1] = v1[2]*v2[0] - v1[0]*v2[2];\
    dest[2] = v1[0]*v2[1] - v1[1]*v2[0];
#define DOT(v1, v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
#define SUB(dest, v1, v2)\
    dest[0]=v1[0]-v2[0];\
    dest[1]=v1[1]-v2[1];\
    dest[2]=v1[2]-v2[2];
*/
inline void CROSS(double dest[3], double v1[3], double v2[3])
{
    dest[0] = v1[1]*v2[2] - v1[2]*v2[1];
    dest[1] = v1[2]*v2[0] - v1[0]*v2[2];
    dest[2] = v1[0]*v2[1] - v1[1]*v2[0];
}

inline double DOT(double v1[3], double v2[3])
{
    return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
}

inline void SUB(double dest[3], double v1[3], double v2[3])
{
    dest[0] = v1[0] - v2[0];
    dest[1] = v1[1] - v2[1];
    dest[2] = v1[2] - v2[2];
}

int intersect_triangle(double orig[3], double dir[3],
        double vert0[3], double vert1[3], double vert2[3],
        double *t, double *u, double *v)
{
    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
    double det, inv_det;

    //find vectors for two edges sharing vert0
    SUB(edge1, vert1, vert0);
    SUB(edge2, vert2, vert0);

    // begin calculating determinant - also used to calculate U parameter
    CROSS(pvec, dir, edge2);

    // if determinant is near zero, ray lies in plane of triangle
    det = DOT(edge1, pvec);

#ifdef TEST_CULL // define TEST_CULL if culling is desired
    if(det < EPSILON)
        return 0;

    // calculate distance from vert0 to ray origin
    SUB(tvec, orig, vert0);

    // calculate U parameter and test bounds
    *u = DOT(tvec, pvec);
    if (*u < 0.0 || *u > det)
        return 0;

    // prepare to test V parameter
    CROSS(qvec, tvec, edge1);

    // calculate V parameter and test bounds
    *v = DOT(dir, qvec);
    if (*v < 0.0 || *u +*v > det)
        return 0;

    // calculate t, scale parameters, ray intersects triangle
    *t = DOT(edge2, qvec);
    int_det = 1.0 / det;
    *t *= inv_det;
    *u *= inv_det;
    *v *= inv_det;
#else           // the non-culling branch
    if (det > -EPSILON && det < EPSILON)
        return 0;
    inv_det = 1.0 / det;

    // calculate distance from vert0 to ray origin
    SUB(tvec, orig, vert0);

    // calculate U parameter and test bounds
    *u = DOT(tvec, pvec) * inv_det;
    if (*u < 0.0 || *u > 1.0)
        return 0;

    // prepare to test V parameters
    CROSS(qvec, tvec, edge1);

    // calculate V paremeter and test bounds
    *v = DOT(dir, qvec) * inv_det;
    if (*v < 0.0 || *u + *v > 1.0)
        return 0;

    //calculate t, ray intersects triangle
    *t = DOT(edge2, qvec) * inv_det;
#endif
    return 1;

}

int main(int argc, char* argv[])
{
    double ori[] = {0.0, 0.0, 0.0};
    double dir[] = {1.0/sqrtf(3), 1.0/sqrtf(3), 1.0/sqrtf(3)};///sqrtf(3)
    double vert0[] = {0.0, 0.0, 1.0};
    double vert1[] = {1.0, 0.0, 0.0};
    double vert2[] = {0.0, 1.0, 0.0};
    double t=0, u, v;

    for (int i=0; i<1; ++i)
    {
        intersect_triangle(ori, dir, vert0, vert1, vert2, &t, &u, &v);
        printf("%lf %lf %lf \n", t, u, v);
    }

    return 0;
}

附文献一篇