Пересечение линии и треугольника в 3D

Question

У меня есть линия и треугольник где-то в трехмерном пространстве. Другими словами, у меня есть 3 точки ([x,y,z] каждая) для треугольника и две точки (также [x,y,z]) для линии.

Мне нужно найти способ, надеюсь, с помощью C++, чтобы выяснить, пересекает ли когда-либо линия треугольник. Линия, параллельная треугольнику и имеющая более одной общей точки, должна рассматриваться как "не пересекающаяся".

Я уже сделал некоторый код, но он не работает, и я всегда получаю ложь, даже когда визуальное представление ясно показывает пересечение.

ofVec3f P1, P2;
P1 = ray.s;
P2 = ray.s + ray.t;

ofVec3f p1, p2, p3;
p1 = face.getVertex(0);
p2 = face.getVertex(1);
p3 = face.getVertex(2);

ofVec3f v1 = p1 - p2;
ofVec3f v2 = p3 - p2;

float a, b, c, d;

a = v1.y * v2.z - v1.z * v2.y;
b = -(v1.x * v2.z - v1.z * v2.x);
c = v1.x * v2.y - v1.y * v2.x;
d = -(a * p1.x + b * p1.y + c * p1.z);

ofVec3f O = P1;
ofVec3f V = P2 - P1;

float t;

t = -(a * O.x + b * O.y + c * O.z + d) / (a * V.x + b * V.y + c * V.z);

ofVec3f p = O + V * t;

float xmin = std::min(P1.x, P2.x);
float ymin = std::min(P1.y, P2.y);
float zmin = std::min(P1.z, P2.z);

float xmax = std::max(P1.x, P2.x);
float ymax = std::max(P1.y, P2.y);
float zmax = std::max(P1.z, P2.z);


if (inside(p, xmin, xmax, ymin, ymax, zmin, zmax)) {
    *result = p.length();
    return true;
}
return false;

А вот и определение inside()

bool primitive3d::inside(ofVec3f p, float xmin, float xmax, float ymin, float ymax, float zmin, float zmax) const {
    if (p.x >= xmin && p.x <= xmax && p.y >= ymin && p.y <= ymax && p.z >= zmin && p.z <= zmax)
        return true;

    return false;
}

Answer 1

1) Если вы просто хотите узнать , пересекает ли линия треугольник (без необходимости пересечения):

Пусть p1,p2,p3 обозначают ваш треугольник

Выберите две точки q1, q2 на линии очень далеко в обоих направлениях.

Пусть SignedVolume(a,b,c,d) обозначает подписанный объем тетраэдра a,b,c,d.

Если SignedVolume(q1,p1,p2,p3) и SignedVolume(q2,p1,p2,p3) имеют разные знаки И SignedVolume(q1,q2,p1,p2), SignedVolume(q1,q2,p2,p3) и SignedVolume(q1,q2,p3,p1) имеют одинаковый знак, то есть пересечение.

SignedVolume(a,b,c,d) = (1/6) * точка (крестик (ba, ca), da)

2) Теперь, если вы хотите пересечение, когда тест в 1) проходит

запишите уравнение прямой в параметрической форме: p(t) = q1 + t*(q2-q1)

Запишите уравнение плоскости: точка (p,N) - точка (p,p1) = 0, где N = крест (p2-p1, p3-p1)

Введите p (t) в уравнение плоскости: точка (q1 + t*(q2-q1), N-p1) = 0

Выведите t = -dot(q1,N-p1)/dot(q1,q2-q1)

Точка пересечения q1 + t * (q2-q1)

Answer 2

@BrunoLevi: ваш алгоритм, похоже, не работает, смотрите следующую реализацию Python:

def intersect_line_triangle(q1,q2,p1,p2,p3):
    def signed_tetra_volume(a,b,c,d):
        return np.sign(np.dot(np.cross(b-a,c-a),d-a)/6.0)

    s1 = signed_tetra_volume(q1,p1,p2,p3)
    s2 = signed_tetra_volume(q2,p1,p2,p3)

    if s1 != s2:
        s3 = signed_tetra_volume(q1,q2,p1,p2)
        s4 = signed_tetra_volume(q1,q2,p2,p3)
        s5 = signed_tetra_volume(q1,q2,p3,p1)
        if s3 == s4 and s4 == s5:
            n = np.cross(p2-p1,p3-p1)
            t = -np.dot(q1,n-p1) / np.dot(q1,q2-q1)
            return q1 + t * (q2-q1)
    return None

Мой тестовый код:

q0 = np.array([0.0,0.0,1.0])
q1 = np.array([0.0,0.0,-1.0])
p0 = np.array([-1.0,-1.0,0.0])
p1 = np.array([1.0,-1.0,0.0])
p2 = np.array([0.0,1.0,0.0])

print(intersect_line_triangle(q0,q1,p0,p1,p2))

дает:

[ 0.  0. -3.] 

вместо ожидаемого

[ 0.  0. 0.]

глядя на линию

t = np.dot(q1,n-p1) / np.dot(q1,q2-q1)

Вычитание p1 из нормали не имеет смысла для меня, вы хотите проецировать от q1 на плоскость треугольника, поэтому вам нужно проецировать вдоль нормали с расстоянием, пропорциональным отношению расстояния от q1 к плоскость и q1-q2 вдоль нормали, верно?

Следующий код исправляет это:

n = np.cross(p2-p1,p3-p1)
t = np.dot(p1-q1,n) / np.dot(q2-q1,n)
return q1 + t * (q2-q1)

Answer 3

Чтобы найти пересечение между линией и треугольником в 3D, используйте следующий подход:

• Вычислить плоскость, поддерживающую треугольник,

• Пересечь линию с плоскостью, поддерживающей треугольник:

  • Если пересечения нет, то пересечения с треугольником нет.

  • Если есть пересечение, убедитесь, что точка пересечения действительно лежит в треугольнике:

    • Каждое ребро треугольника вместе с нормалью плоскости, поддерживающей треугольник, определяет полупространство, ограничивающее внутреннюю часть треугольника (соответствующая ограничивающая плоскость может быть получена из вершин нормали и ребер),
    • Убедитесь, что точка пересечения лежит внутри всех краевых полупространств.

Вот пример кода с подробными вычислениями, которые должны работать:

// Compute the plane supporting the triangle (p1, p2, p3)
//     normal: n
//     offset: d
//
// A point P lies on the supporting plane iff n.dot(P) + d = 0
//
ofVec3f v21 = p2 - p1;
ofVec3f v31 = p3 - p1;

ofVec3f n = v21.getCrossed(v31);
float d = -n.dot(p1);

// A point P belongs to the line from P1 to P2 iff
//     P = P1 + t * (P2 - P1)
//
// Find the intersection point P(t) between the line and
// the plane supporting the triangle:
//     n.dot(P) + d = 0
//                  = n.dot(P1 + t (P2 - P1)) + d
//                  = n.dot(P1) + t n.dot(P2 - P1) + d
//
//     t = -(n.dot(P1) + d) / n.dot(P2 - P1)
//
ofVec3f P21 = P2 - P1;
float nDotP21 = n.dot(P21);

// Ignore line parallel to (or lying in) the plane
if (fabs(nDotP21) < Epsilon)
    return false;

float t = -(n.dot(P1) + d) / nDotP21;
ofVec3f P = P1 + t * P21;

// Plane bounding the inside half-space of edge (p1, p2): 
//     normal: n21 = n x (p2 - p1)
//     offset: d21 = -n21.dot(p1)
//
// A point P is in the inside half-space iff n21.dot(P) + d21 > 0
//

// Edge (p1, p2)
ofVec3f n21 = n.cross(v21);
float d21 = -n21.dot(p1);

if (n21.dot(P) + d21 <= 0)
    return false;

// Edge (p2, p3)
ofVec3f v32 = p3 - p2;
ofVec3f n32 = n.cross(v32);
float d32 = -n32.dot(p2);

if (n32.dot(P) + d32 <= 0)
    return false;

// Edge (p3, p1)
ofVec3f n13 = n.cross(-v31);
float d13 = -n13.dot(p3);

if (n13.dot(P) + d13 <= 0)
    return false;

return true;

Некоторые комментарии к коду размещены с вопросом:

• Предопределенные операции ofVec3f (.dot() а также .cross() для геометрических продуктов и т. д.) следует отдавать предпочтение при их наличии (более читабельно, избегать ошибок при реализации и т. д.),
• Код вначале следует описанному выше подходу, но затем проверяет только то, что точка пересечения находится в трехмерной ограничивающей рамке отрезка отрезка [P1, P2]. Это в сочетании с возможными другими ошибками может объяснить, почему результаты неверны.
• Можно убедиться, что точка пересечения находится в трехмерном ограничивающем прямоугольнике (целого) треугольника. Хотя этого недостаточно, чтобы гарантировать пересечение, его можно использовать для того, чтобы отбирать точки, которые явно не пересекаются, и избегать дальнейших сложных вычислений.

Answer 4

У меня есть другой способ сделать это, который, как я обнаружил в моем рендерере, намного быстрее, чем первый способ, предложенный БруноЛеви. (Я не реализовал второй способ)

Точки A, B, C — вершины треугольника
O — начало луча
D — направление луча (не требует нормализации, просто ближе к началу координат, чем треугольник)

Проверить, находится ли направление (D+O) внутри тетраэдра A, B, C, O

      bool SameSide(vec3 A, vec3 B, vec3 C, vec3 D, vec3 p)
{
    vec3 normal = cross(B - A, C - A);
    float dotD = dot(normal, D - A);
    float dotP = dot(normal, p - A);
    return signbit(dotD) == signbit(dotP);
}

bool LineIntersectTri(vec3 A, vec3 B, vec3 C, vec3 O, vec3 D)
{
    return SameSide(A, B, C, O, O+D) &&
           SameSide(B, C, O, A, O+D) &&
           SameSide(C, O, A, B, O+D) &&
           SameSide(O, A, B, C, O+D);               
}

Если D меняется, а все остальное остается прежним (например, в рендерере raycasting), тогда нормальный и dotP не нужно пересчитывать; Вот почему я нашел это намного быстрее

Код пришел из этого ответа /questions/21231818/kak-proverit-nahoditsya-li-tochka-v-tetraedre-ili-net/21231831#21231831