c - 算法错误地说射线与其上方的三角形相交

Question

这是许多类似的射线-三角形相交算法之一。我测试过的所有其他算法也对这些数字返回 true，而光线显然没有穿过三角形。光线从 y=0 到 y=1，而三角形在 y = 2.3 上是平坦的。

这不是一个缠绕问题，因为它永远不会返回真(缠绕问题会解释漏报，而不是误报)。

在 C 或 C++ 中重现所需的所有代码都包含在此处。

我错过了什么？

#define vector(a,b,c) \
(a)[0] = (b)[0] - (c)[0];   \
(a)[1] = (b)[1] - (c)[1];   \
(a)[2] = (b)[2] - (c)[2];

#define crossProduct(a,b,c) \
(a)[0] = (b)[1] * (c)[2] - (c)[1] * (b)[2]; \
(a)[1] = (b)[2] * (c)[0] - (c)[2] * (b)[0]; \
(a)[2] = (b)[0] * (c)[1] - (c)[0] * (b)[1];

#define innerProduct(v,q) \
((v)[0] * (q)[0] + \
(v)[1] * (q)[1] + \
(v)[2] * (q)[2])

#define DOT(A,B) \
((A)[0] * (B)[0] + (A)[1] * (B)[1] + (A)[2] * (B)[2])


    int intersect3D_RayTriangle( )
{
    //    dir, w0, w;           // ray vectors
    double     r, a, b;              // params to calc ray-plane intersect

    // output: Point* I

    //Ray R
    double origin[3] = {0,0,0};//{orig[0],orig[1],orig[2]};
    double direction[3] = {0,1,0};//{dir[0],dir[1],dir[2]};

    //Triangle T
    double corner1[3] = {3,    2.3,    -4 };//{v0[0],v0[1],v0[2]};
    double corner2[3] = {-7,   2.3,    2};//{v1[0],v1[1],v1[2]};
    double corner3[3] = {3,    2.3,    2};// v2[0],v2[1],v2[2]};

    //     Vector    u, v, n;              // triangle vectors
    double u[3] = {corner2[0]-corner1[0],corner2[1]-corner1[1],corner2[2]-corner1[2]};
    double v[3] = {corner3[0]-corner1[0],corner3[1]-corner1[1],corner3[2]-corner1[2]};
    double n[3] = {0,0,0};
    double e1[3],e2[3],h[3],q[3];
    double f;



    // get triangle edge vectors and plane normal
    crossProduct(n, u, v);
    if ((n[0] == 0) && (n[1] == 0) && (n[2] == 0))             // triangle is wonky
        return -1;                  // do not deal with this case

    // dir = R.P1 - R.P0;              // ray direction vector
    double rayDirection[3] = {direction[0] - origin[0], direction[1] - origin[1], direction[2] - origin[2]};

    //w0 = R.P0 - T.V0;
    double w0[3]  = {origin[0] - corner1[0], origin[1] - corner1[1], origin[2] - corner1[2]};

    a = -DOT(n,w0);
    b = DOT(n,rayDirection);
    if (fabs(b) < __DBL_EPSILON__) {     // ray is  parallel to triangle plane
        if (a == 0)                 // ray lies in triangle plane
            return 2;
        else return 0;              // ray disjoint from plane
    }

    // get intersect point of ray with triangle plane
    r = a / b;
    if (r < 0.0)                    // ray goes away from triangle
        return 0;                   // => no intersect
    // for a segment, also test if (r > 1.0) => no intersect

    //*I = R.P0 + r * dir;            // intersect point of ray and plane
    double I[3] = {0,0,0};
    I[0] = origin[0] + rayDirection[0] * r;
    I[1] = origin[1] + rayDirection[1] * r;
    I[2] = origin[2] + rayDirection[2] * r;

    // is I inside T?
    double    uu, uv, vv, wu, wv, D;
    uu = DOT(u,u);
    uv = DOT(u,v);
    vv = DOT(v,v);

    double w[3] = {0,0,0};
    w[0] = I[0] - corner1[0];
    w[1] = I[1] - corner1[1];
    w[2] = I[2] - corner1[2];


    wu = DOT(w,u);
    wv = DOT(w,v);
    D = uv * uv - uu * vv;

    // get and test parametric coords
    double s, t;
    s =  (uv * wv - vv * wu) / D;
    if (s < 0.0 || s > 1.0)         // I is outside T
        return 0;
    t = (uv * wu - uu * wv) / D;
    if (t < 0.0 || (s + t) > 1.0)  // I is outside T
        return 0;

    return 1;                       // I is in T
}

Answer 1

代码适用于“射线”。

OP 期望“射线”代码的功能类似于“段”代码。

可以使用 r 值来测试“段”排除。

if (r > 1.0) return 0;