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

问题描述 投票：0回答：5

我有一个由 3d 空间中的 3 个点定义的三角形。我还有一个由 3d 空间中的 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 交集算法。请参阅下面的替代解决方案的答案。

c++ math 3d polygon
5个回答
5
投票

这是解决您问题的一种方法。计算四面体的体积 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) 相交。

如何根据角坐标计算四面体的体积 web，以及C中的计算几何。

3
投票

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

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

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

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

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)
2
投票

感谢您的帮助！这是一个替代解决方案。问题是针对 c++ 的，正如 meowgoestedog 指出的那样，我可以尝试使用 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;
}
0
投票

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;
      }
    }
  }
0
投票

Unity C# 版本：

public static bool RaycastTriangle(
    Vector3 origin, 
    Vector3 pointB, 
    Vector3 a, 
    Vector3 b, 
    Vector3 c, 
    out Vector3 intersection
) {
    intersection = Vector3.zero;
    var direction = pointB - origin;
    var normalizedDirection = direction.normalized;
    
    Vector3 edge1 = b - a;
    Vector3 edge2 = c - a;
    Vector3 normal = Vector3.Cross(normalizedDirection, edge2);

    float dotProduct = Vector3.Dot(edge1, normal);

    if (dotProduct > -float.Epsilon && dotProduct < float.Epsilon)
        return false; // Parallel

    float inverseDotProduct = 1f / dotProduct;
    Vector3 toStart = origin - a;
    float triangleParam = inverseDotProduct * Vector3.Dot(toStart, normal);

    if (triangleParam is < 0f or > 1f) return false;
    Vector3 edgeCross = Vector3.Cross(toStart, edge1);
    float raycastParam = inverseDotProduct * Vector3.Dot(normalizedDirection, edgeCross);

    if (raycastParam < 0f || triangleParam + raycastParam > 1f) return false;
    
    float distance = inverseDotProduct * Vector3.Dot(edge2, edgeCross);
    if (distance > float.Epsilon && distance < direction.magnitude)
    {
        intersection = origin + normalizedDirection * distance;
        return true;
    }
    return false;
}
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.