TriangleMesh - is_self_intersecting return wrong result

I do see similar problem that not all intersecting triangles are reported. When removing the intersecting triangles and removing all small clusters of triangles not connected any more from the surface does not clear out all intersections. especially in very narrow surface loops it can happen that intersecting triangles remain.

One reason for this is the test in lines 1360 to 1367 of TriangleMesh.cpp in function TriangleMesh::GetSelfIntersectingTriangles

  // check if neighbour triangle
            if (tria_p(0) == tria_q(0) || tria_p(0) == tria_q(1) ||
                tria_p(0) == tria_q(2) || tria_p(1) == tria_q(0) ||
                tria_p(1) == tria_q(1) || tria_p(1) == tria_q(2) ||
                tria_p(2) == tria_q(0) || tria_p(2) == tria_q(1) ||
                tria_p(2) == tria_q(2)) {
                continue;
            }

This test rules out any pair of triangles which share at least one point instead of ruling out triangles sharing a common edge only. As two triangles which share only one point without having any edge in common can easily intersect each other. They just need to share edges with an intermediate triangle having the same point in common.

Concerning the issues the others having: As far as I do understand the code of the IntersectionTest::TriangleTriangle3d method in IntersectionTest.cpp file

bool IntersectionTest::TriangleTriangle3d(const Eigen::Vector3d& p0,
                                          const Eigen::Vector3d& p1,
                                          const Eigen::Vector3d& p2,
                                          const Eigen::Vector3d& q0,
                                          const Eigen::Vector3d& q1,
                                          const Eigen::Vector3d& q2) {
    const Eigen::Vector3d mu = (p0 + p1 + p2 + q0 + q1 + q2) / 6;
    const Eigen::Vector3d sigma =
            (((p0 - mu).array().square() + (p1 - mu).array().square() +
              (p2 - mu).array().square() + (q0 - mu).array().square() +
              (q1 - mu).array().square() + (q2 - mu).array().square()) /
             5)
                    .sqrt() +
            1e-12;
    Eigen::Vector3d p0m = (p0 - mu).array() / sigma.array();
    Eigen::Vector3d p1m = (p1 - mu).array() / sigma.array();
    Eigen::Vector3d p2m = (p2 - mu).array() / sigma.array();
    Eigen::Vector3d q0m = (q0 - mu).array() / sigma.array();
    Eigen::Vector3d q1m = (q1 - mu).array() / sigma.array();
    Eigen::Vector3d q2m = (q2 - mu).array() / sigma.array();
    return NoDivTriTriIsect(p0m.data(), p1m.data(), p2m.data(), q0m.data(),
                            q1m.data(), q2m.data()) != 0;
}

It standardizes both triangles to mean vertex (0,0,0) and standardises the length of edge vectors to the standard deviation of the edge lenghs if all edge vectors and than assumingly applies some machine learning to guess whether the two triangles intersect or not or other not sufficiently deterministic approach to do so. Given the multiple flaws of your current approach to triangle and thus Triangle mesh self intersection i do suggest you resort to classical edge plane intersection and barycentric on triangle testing or replace current approach by any othe robust and reliable intersection test.

// these yield penetration distance independent whether edge runs in direction of triangle normal or oposit to it.
k1 = dot(p1 - q1, n1) / dot( ( q2-q1)/|q2 - q1|,n1)  
k2 = dot(p1 - q2, n1) / dot( (q3 - q2)/|q3-q2|,n1)
k3 = dot(p1 - q3, n1) / dot( (q1-q3)/|q1-q3|,n1)
l1 = dot(q1 - p1, n2) / dot( ( p2-p1)/|p2 - p1|,n2) 
l2 = dot(q1 - p2, n2) / dot( (p3 - p2)/|p3-p2|,n2)
l3 = dot(q1 - p3, n2) / dot( (p1-p3)/|p1-p3|,n2)

// test if the penetrations distances are within the edge length of the edge tested which is the case if penetrating independent of edge orientation with respect to triangle normal
if ( ( 0<= k1 <= |q2-q1|) or (0<= k2 <= |q3-q2|) or ( 0<=k3<=|q1-q3|)  or (0<= l1 <= |p2-p1|) or (0<= l2 <= |p3-p2|) or ( 0<=l3<=|p1-p3|) ) {
   (q1',q2',q3',p1',p2',p3') = (q1,q2,q3,p1,p2,p3) + (k1,k2,k3,l1,l2,l3) * <edgevectors>
   is_on_tri(q1',q2',q3' ,p1,p2,p3) , is_on_tri(p1',p2',p3',q1,q2,q3)

and the barycentric check if point is on triangle is

is_on_tri(qX',p1,p2,p3)
A = |cross(p2-p1,p3-p1)| // should be memorized to avoid multiple recomputations
alpha=|cross(p2-qX',p3-qX')| / A
if 0<= alpha <= 1:
    beta = |cross(p3-qX',p1-qX')|/A
    if 0<=beta<=1 and 0<= (alpha+beta= <= 1): // note testing if alpha + beta is larger than 1 alone is not sufficient
         gamma = |cross(p1-qX',p2-qX')|/A         // the areas of the first two sub triangles could be together smaller than A
         return 0<= gamma <= 1 and 0<=gamma+alpha+beta<=1 // even though qX' is outside the triangle
 return false

This might not be the fastest amongst available robust intersection tests but can cope with all kinds of triangle triangle intersections including those your approach is failing.

Reading the paper where the author or the test explains his approach than i do miss the normalization of the plane normal to 1 for properly computing the intersection of edge and triangle plane using plane equation and standardization of points or edges does not yield the situation that both vectors have length 1 it just ensures that when computing standard deviation of the normalized edge lenght it turns out to be 1. The individual edges still may have edge length different from 1. And shifting origin of all vertices to their mean doesn't change this either. And in his paper the author does only talk about special cases but does not give any indication how he tests if plane intersection point lies within triangle 1. In other words i do doubt that his optimizations are capable to capture all possible cases, i guess several are missed or points falling outside the triangle as observed by OP and others are considered on triangle.

EDIT Reading the Triangle Ray intersection paper there is the proof and reason why the OP and others observe the strange behavior. Moller uses Barycentric coordinates. And for his equation 2 he implicitly takes the assumption that it is already proven that the ray intersects the triangle plane at a point within the triangle. therefore he only computes the Barycentric coordinates u and v and uses the Barycentric condition s + u + v = 1 to derive s from u and v. In general case where ray triangle intersection test should be used to test if point where ray intersects with triangle plane is located inside triangle this is insufficient. There are cases where 0<=u<=1, 0<=v<=1 and 0<=u+v<=1 but the intersection point is located outside the triangle. And these cases the equation two Akenine Moller uses does not capture properly as it sit still possible that 0>s or s>1 or s+u+v>1 when the other four conditions are met. Only when all conditions are met than the intersection point will be a point on the triangle.

1) 0<=u<=1 // not on triangle if fails
2) 0<=v<=1 // not on triangle if fails
3) 0<=u+v<=1 // not on triangle if fails
4) 0<=s<=1 // not on triangle if false
5) 0<=s+u+v<=1 // on triangle only if this and all the afore listed succeed

The order is by will as v only needs to be computed if 1) succeeds and s only needs to be computed 1) 2) and 3) succeed.

To be hones i also made this mistake when starting coding edge triangle intersections manually observing strange behaviors like the OP until until accepting.

isl-org / Open3D

TriangleMesh - is_self_intersecting return wrong result #5117

Checklist

Describe the issue

Steps to reproduce the bug

Error message

Expected behavior

Open3D, Python and System information

Additional information