2017-07-10 2 views
2

これは、多くの同様のレイ・トライアングル交差アルゴリズムの1つです。私がテストした他のアルゴリズムもすべて、これらの数値に対して真を返しますが、レイは明らかに三角形を横切っていません。光線は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 
} 
+0

こんにちはチルトン、ようこそ! 「良い」質問をするには、次のリンクを必ず読んでください。それは、人々が関連性のある回答でよりうまく答えるのを助けるでしょう。 https://stackoverflow.com/help/how-to-ask –

+0

Noted、edited。ありがとうございました。 –

+0

好奇心。なぜ主に「ダブル」を使用するのですが、ドットプロダクトは「浮動」ですか? – chux

答えて

1

"rays"のコードは正常に機能します。

"光線"コードは "セグメント"のように機能すると予想されました。

r値を「セグメント」除外のテストに使用できます。

if (r > 1.0) return 0; 
+0

ありがとう、それは完璧に動作します! –

関連する問題