如何知道线段是否与 3d space 中的三角形相交？

Question

我有一个由 3d 中的 3 个点定义的三角形 space。我还有一个由 3d space 中的 2 个点定义的线段。我想知道它们是否相交。我真的不需要知道交点。

我不懂微积分，但我懂一些三角函数。我知道一些关于矩阵的知识，但我很了解向量（特别是 3d 向量）。请保持简单。

你能帮我完成示例问题吗：

三角形：

a: -4, 3, 0

b: 4, 3, 0

c: -3, -5, 4

线段：

d: 1, -2, 0

e: -2, 6, 2

编辑：

我打算在 C++ 物理引擎中使用它。

一个答案涉及从 4 个顶点计算四面体体积。请提供公式或在代码中显示。

更新：

正如 meowgoesthedog 指出的，我可以尝试使用 Moller-Trumbore 交集算法。请参阅下面的替代解决方案的答案。

Answer 1

这是解决您问题的一种方法。计算四面体的体积 Td = (a,b,c,d) 和 Te = (a,b,c,e)。如果 Td 或 Te 的体积为零，则线段 de 位于包含三角形 (a,b,c) 的平面上。如果 Td 和 Te 的体积具有相同的符号，那么de严格偏向一侧，没有交集。如果 Td 和 Te 取反符号，然后 de 穿过包含 (a,b,c) 的平面。

从那里有几种策略。一种是计算 de 交叉点 p 那架飞机然后投影到2D，解决2D中的三角点问题

另一种方法是计算四面体 (a,b,d,e)、(b,c,d,e) 和 (c,a,d,e) 的体积。那么只有当三个符号都相同时，才与三角形 (a,b,c) 相交。

如何根据角坐标计算四面体的体积网络，也在 Computational Geometry in C.

Answer 2

我实施了约瑟夫在 python 中给出的出色答案，并认为我会分享。该函数采用一组线段和三角形，并计算每条线段是否与任何给定三角形相交。

函数的第一个输入是一个 2xSx3 线段数组，其中第一个索引指定线段的起点或终点，第二个索引指的是第 s^ 个线段，第三个索引指向线段点的 x、y、z 坐标。

第二个输入是三角形顶点的 3XTX3 数组，其中第一个索引指定三个顶点之一（不必按任何特定顺序），第二个索引指的是第 t^th 个三角形, 第三个索引指向三角形顶点的 x,y,z 坐标。

此函数的输出是一个大小为 S 的二进制数组，它告诉您第 s^ 条线段是否与给定的任何三角形相交。如果你想知道线段与哪些三角形相交，那么只需删除函数最后一行的总和即可。

def signedVolume(a, b, c, d):
    """Computes the signed volume of a series of tetrahedrons defined by the vertices in 
    a, b c and d. The ouput is an SxT array which gives the signed volume of the tetrahedron defined
    by the line segment 's' and two vertices of the triangle 't'."""

    return np.sum((a-d)*np.cross(b-d, c-d), axis=2)

def segmentsIntersectTriangles(s, t):
    """For each line segment in 's', this function computes whether it intersects any of the triangles
    given in 't'."""
    # compute the normals to each triangle
    normals = np.cross(t[2]-t[0], t[2]-t[1])
    normals /= np.linalg.norm(normals, axis=1)[:, np.newaxis]

    # get sign of each segment endpoint, if the sign changes then we know this segment crosses the
    # plane which contains a triangle. If the value is zero the endpoint of the segment lies on the 
    # plane.
    # s[i][:, np.newaxis] - t[j] -> S x T x 3 array
    sign1 = np.sign(np.sum(normals*(s[0][:, np.newaxis] - t[2]), axis=2)) # S x T
    sign2 = np.sign(np.sum(normals*(s[1][:, np.newaxis] - t[2]), axis=2)) # S x T

    # determine segments which cross the plane of a triangle. 1 if the sign of the end points of s is 
    # different AND one of end points of s is not a vertex of t
    cross = (sign1 != sign2)*(sign1 != 0)*(sign2 != 0) # S x T 

    # get signed volumes
    v1 = np.sign(signedVolume(t[0], t[1], s[0][:, np.newaxis], s[1][:, np.newaxis])) # S x T
    v2 = np.sign(signedVolume(t[1], t[2], s[0][:, np.newaxis], s[1][:, np.newaxis])) # S x T
    v3 = np.sign(signedVolume(t[2], t[0], s[0][:, np.newaxis], s[1][:, np.newaxis])) # S x T

    same_volume = np.logical_and((v1 == v2), (v2 == v3)) # 1 if s and t have same sign in v1, v2 and v3

    return (np.sum(cross*same_volume, axis=1) > 0)

Answer 3

感谢您的帮助！这是一个替代解决方案。问题是针对 c++ 的，正如 meowgoesthedog 指出的那样，我可以尝试使用 Moller-Trumbore 交集算法。这是我想出的：

#include <math.h>

class vec3 {
public:
    float x, y, z;

    float dot(const vec3 & b) {
        return vec3::x * b.x + vec3::y * b.y + vec3::z * b.z;
    }

    vec3 cross(const vec3 & b) {
        return vec3::vec3(
            vec3::y * b.z - vec3::z * b.y,
            vec3::z * b.x - vec3::x * b.z,
            vec3::x * b.y - vec3::y * b.x
        );
    }

    vec3 normalize() {
        const float s = 1.0f / sqrtf(vec3::x * vec3::x + vec3::y * vec3::y + vec3::z * vec3::z);
        return vec3::vec3(vec3::x * s, vec3::y * s, vec3::z * s);
    }

    vec3 operator+(const vec3 & b) {
        return vec3::vec3(
            vec3::x + b.x,
            vec3::y + b.y,
            vec3::z + b.z
        );
    }
    vec3 operator+=(const vec3 & b) {
        *this = vec3::operator+(b);
        return *this;
    }

    vec3 operator-(const vec3 & b) {
        return vec3::vec3(
            vec3::x - b.x,
            vec3::y - b.y,
            vec3::z - b.z
        );
    }
    vec3 operator-=(const vec3 & b) {
        *this = vec3::operator-(b);
        return *this;
    }

    vec3 operator*(const vec3 & b) {
        return vec3::vec3(
            vec3::x * b.x,
            vec3::y * b.y,
            vec3::z * b.z
        );
    }
    vec3 operator*=(const vec3 & b) {
        *this = vec3::operator*(b);
        return *this;
    }
    vec3 operator*(float b) {
        return vec3::vec3(
            vec3::x * b,
            vec3::y * b,
            vec3::z * b
        );
    }
    vec3 operator*=(float b) {
        *this = vec3::operator*(b);
        return *this;
    }

    vec3 operator/(const vec3 & b) {
        return vec3::vec3(
            vec3::x / b.x,
            vec3::y / b.y,
            vec3::z / b.z
        );
    }
    vec3 operator/=(const vec3 & b) {
        *this = vec3::operator/(b);
        return *this;
    }
    vec3 operator/(float b) {
        return vec3::vec3(
            vec3::x * b,
            vec3::y * b,
            vec3::z * b
        );
    }
    vec3 operator/=(float b) {
        *this = vec3::operator/(b);
        return *this;
    }

    vec3(float x, float y, float z) {
        vec3::x = x;
        vec3::y = y;
        vec3::z = z;
    }
    vec3(float x) {
        vec3::x = x;
        vec3::y = x;
        vec3::z = x;
    }
    vec3() {
        //
    }
    ~vec3() {
        //
    }
};

#define EPSILON 0.000001f
bool lineSegIntersectTri(
    vec3 line[2],
    vec3 tri[3],
    vec3 * point
) {
    vec3 e0 = tri[1] - tri[0];
    vec3 e1 = tri[2] - tri[0];

    vec3 dir = line[1] - line[0];
    vec3 dir_norm = dir.normalize();

    vec3 h = dir_norm.cross(e1);
    const float a = e0.dot(h);

    if (a > -EPSILON && a < EPSILON) {
        return false;
    }

    vec3 s = line[0] - tri[0];
    const float f = 1.0f / a;
    const float u = f * s.dot(h);

    if (u < 0.0f || u > 1.0f) {
        return false;
    }

    vec3 q = s.cross(e0);
    const float v = f * dir_norm.dot(q);

    if (v < 0.0f || u + v > 1.0f) {
        return false;
    }

    const float t = f * e1.dot(q);

    if (t > EPSILON && t < sqrtf(dir.dot(dir))) { // segment intersection
        if (point) {
            *point = line[0] + dir_norm * t;
        }

        return true;
    }

    return false;
}

对于运行一些测试：

#include <stdio.h>

const char * boolStr(bool b) {
    if (b) {
        return "true";
    }

    return "false";
}

int main() {
    vec3 tri[3] = {
        { -1.0f, -1.0f, 0.0f },
        { 1.0f, -1.0f, 0.0f },
        { 1.0f, 1.0f, 0.0f },
    };

    vec3 line0[2] = { // should intersect
        { 0.5f, -0.5f, -1.0f },
        { 0.5f, -0.5f, 1.0f },
    };

    vec3 line1[2] = { // should not intersect
        { -0.5f, 0.5f, -1.0f },
        { -0.5f, 0.5f, 1.0f },
    };

    printf(
        "line0 intersects? : %s\r\n"
        "line1 intersects? : %s\r\n",
        boolStr(lineSegIntersectTri(line0, tri, NULL)),
        boolStr(lineSegIntersectTri(line1, tri, NULL))
    );

    return 0;
}

Answer 4

C# 版本：

public class AlgoritmoMollerTrumbore
  {
    private const double EPSILON = 0.0000001;

    public static bool lineIntersectTriangle(Point3D[] line,
                                             Point3D[] triangle,
                                             out Point3D outIntersectionPoint)
    {
      outIntersectionPoint = new Point3D(0, 0, 0);

      Point3D rayOrigin = line[0];
      Vector3D rayVector = Point3D.Subtract(line[1], line[0]);
      rayVector.Normalize();

      Point3D vertex0 = triangle[0];
      Point3D vertex1 = triangle[1]; 
      Point3D vertex2 = triangle[2];

      Vector3D edge1 = Point3D.Subtract(vertex1, vertex0);
      Vector3D edge2 = Point3D.Subtract(vertex2, vertex0);
      Vector3D h = Vector3D.CrossProduct(rayVector, edge2);

      double a = Vector3D.DotProduct(edge1, h);
      if (a > -EPSILON && a < EPSILON)
      {
        return false;    // This ray is parallel to this triangle.
      }

      double f = 1.0 / a;      
      Vector3D s = Point3D.Subtract(rayOrigin, vertex0);

      double u = f * (Vector3D.DotProduct(s, h));
      if (u < 0.0 || u > 1.0)
      {
        return false;
      }
      Vector3D q = Vector3D.CrossProduct(s, edge1);
      double v = f * Vector3D.DotProduct(rayVector, q);
      if (v < 0.0 || u + v > 1.0)
      {
        return false;
      }
      // At this stage we can compute t to find out where the intersection point is on the line.
      double t = f * Vector3D.DotProduct(edge2, q);
      if (t > EPSILON && t < Math.Sqrt(Vector3D.DotProduct(rayVector, rayVector))) // ray intersection
      {
        outIntersectionPoint = rayOrigin + rayVector * t;
        return true;
      }
      else // This means that there is a line intersection but not a ray intersection.
      {
        return false;
      }
    }
  }

如何知道线段是否与 3d space 中的三角形相交？

How to know IF a line segment intersects a triangle in 3d space?

c++

math

3d

polygon