Fix false positive triangle-triangle intersections

CHANGELOG.md

@@ -1,4 +1,5 @@
 ## Main
+-   Fix rotation-dependent false-positive triangle-triangle intersection results (PR #7539) (issue #5117).
 -   Add symmetric ICP registration to the legacy and Tensor pipelines (PR #7276).
 -   Replace OpenMP with oneAPI TBB for all CPU parallelism; Open3D no longer depends on OpenMP. This removes the `libomp` / `libgomp` runtime dependency and the thread oversubscription and crashes caused by loading multiple OpenMP runtimes in one process (e.g. alongside PyTorch in Python). The `WITH_OPENMP` CMake option is removed, oneTBB >= 2021.4.0 is required, and `OMP_NUM_THREADS` is replaced by `open3d.utility.set_max_threads()` (C++: `utility::SetMaxThreads()` or a `tbb::task_arena`). `utility::OMPProgressBar` is removed in favor of the thread-safe `utility::ProgressBar`; `utility::GetThreadNum()` and `utility::InParallel()` are removed (PR #6626) (issues #6196, #6544, #6750)
 -   Add point cloud smoothing algorithms: Moving Least Squares (MLS), Laplacian, Taubin, and bilateral smoothing. These methods provide flexible noise reduction for point clouds with different preservation characteristics (PR #7419).

cpp/open3d/geometry/IntersectionTest.cpp

@@ -14,6 +14,39 @@
 
 namespace open3d {
 namespace geometry {
+namespace {
+
+constexpr double kUnitNormalPlaneDistanceTolerance = 1e-6;
+
+// A triangle cannot intersect another triangle when all of its vertices are
+// unambiguously on the same side of the other triangle's supporting plane.
+bool SeparatedByTrianglePlane(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 normal = (p1 - p0).cross(p2 - p0);
+    const double normal_norm = normal.norm();
+    if (!std::isfinite(normal_norm) || normal_norm == 0.0) {
+        return false;
+    }
+    const Eigen::Vector3d unit_normal = normal / normal_norm;
+    const double d0 = unit_normal.dot(q0 - p0);
+    const double d1 = unit_normal.dot(q1 - p0);
+    const double d2 = unit_normal.dot(q2 - p0);
+    if (!std::isfinite(d0) || !std::isfinite(d1) || !std::isfinite(d2)) {
+        return false;
+    }
+    return (d0 > kUnitNormalPlaneDistanceTolerance &&
+            d1 > kUnitNormalPlaneDistanceTolerance &&
+            d2 > kUnitNormalPlaneDistanceTolerance) ||
+           (d0 < -kUnitNormalPlaneDistanceTolerance &&
+            d1 < -kUnitNormalPlaneDistanceTolerance &&
+            d2 < -kUnitNormalPlaneDistanceTolerance);
+}
+
+}  // namespace
 
 bool IntersectionTest::AABBAABB(const Eigen::Vector3d& min0,
                                 const Eigen::Vector3d& max0,
@@ -51,6 +84,10 @@ bool IntersectionTest::TriangleTriangle3d(const Eigen::Vector3d& p0,
     Eigen::Vector3d q0m = (q0 - mu).array() / sigma.array();
     Eigen::Vector3d q1m = (q1 - mu).array() / sigma.array();
     Eigen::Vector3d q2m = (q2 - mu).array() / sigma.array();
+    if (SeparatedByTrianglePlane(p0m, p1m, p2m, q0m, q1m, q2m) ||
+        SeparatedByTrianglePlane(q0m, q1m, q2m, p0m, p1m, p2m)) {
+        return false;
+    }
     return NoDivTriTriIsect(p0m.data(), p1m.data(), p2m.data(), q0m.data(),
                             q1m.data(), q2m.data()) != 0;
 }

cpp/tests/geometry/IntersectionTest.cpp

@@ -7,11 +7,75 @@
 
 #include "open3d/geometry/IntersectionTest.h"
 
+#include <array>
+#include <cmath>
+
 #include "tests/Tests.h"
 
 namespace open3d {
 namespace tests {
 
+TEST(IntersectionTest, TriangleTriangle3dIssue5117) {
+    const std::array<Eigen::Vector3d, 3> p{
+            Eigen::Vector3d(0.0, 0.13918686, 1.0),
+            Eigen::Vector3d(0.0, 0.0, 1.1270161),
+            Eigen::Vector3d(1.0, 0.0, 1.0284119)};
+    const std::array<Eigen::Vector3d, 3> q{
+            Eigen::Vector3d(1.0, 1.1269569, 0.0),
+            Eigen::Vector3d(1.0, 0.03113556, 1.0),
+            Eigen::Vector3d(2.0, 1.0189056, 0.0)};
+
+    // These separated triangles used to be reported as intersecting until a
+    // rigid rotation moved their plane values beyond an absolute tolerance.
+    const auto expect_separated = [](const auto& p_vertices,
+                                     const auto& q_vertices) {
+        for (int p_offset = 0; p_offset < 3; ++p_offset) {
+            for (int q_offset = 0; q_offset < 3; ++q_offset) {
+                EXPECT_FALSE(geometry::IntersectionTest::TriangleTriangle3d(
+                        p_vertices[p_offset], p_vertices[(p_offset + 1) % 3],
+                        p_vertices[(p_offset + 2) % 3], q_vertices[q_offset],
+                        q_vertices[(q_offset + 1) % 3],
+                        q_vertices[(q_offset + 2) % 3]));
+                EXPECT_FALSE(geometry::IntersectionTest::TriangleTriangle3d(
+                        q_vertices[q_offset], q_vertices[(q_offset + 1) % 3],
+                        q_vertices[(q_offset + 2) % 3], p_vertices[p_offset],
+                        p_vertices[(p_offset + 1) % 3],
+                        p_vertices[(p_offset + 2) % 3]));
+            }
+        }
+    };
+    expect_separated(p, q);
+
+    const double sqrt_half = std::sqrt(0.5);
+    Eigen::Matrix3d rotation;
+    rotation << 1.0, 0.0, 0.0, 0.0, sqrt_half, sqrt_half, 0.0, -sqrt_half,
+            sqrt_half;
+    for (double scale : {1e-3, 1.0, 1e3}) {
+        std::array<Eigen::Vector3d, 3> transformed_p;
+        std::array<Eigen::Vector3d, 3> transformed_q;
+        const Eigen::Vector3d translation(10.0, -20.0, 30.0);
+        for (int i = 0; i < 3; ++i) {
+            transformed_p[i] = scale * rotation * p[i] + translation;
+            transformed_q[i] = scale * rotation * q[i] + translation;
+        }
+        expect_separated(transformed_p, transformed_q);
+    }
+}
+
+TEST(IntersectionTest, TriangleTriangle3dNearParallelIntersection) {
+    // Keep uncertain near-parallel cases in the legacy predicate: these two
+    // triangles cross even though their normals are almost parallel.
+    const Eigen::Vector3d p0(-1.0, -1.0, -2.0);
+    const Eigen::Vector3d p1(1.0, -1.0, 0.0);
+    const Eigen::Vector3d p2(0.0, 1.0, 1.0);
+    const Eigen::Vector3d q0(-1.0, -1.0, -2.000002);
+    const Eigen::Vector3d q1(1.0, -1.0, 0.000002);
+    const Eigen::Vector3d q2(0.0, 1.0, 1.0);
+
+    EXPECT_TRUE(geometry::IntersectionTest::TriangleTriangle3d(p0, p1, p2, q0,
+                                                               q1, q2));
+}
+
 TEST(IntersectionTest, PointsCoplanar) {
     Eigen::Vector3d p0(0, 0, 0);
     Eigen::Vector3d p1(1, 0, 0);

cpp/tests/geometry/TriangleMesh.cpp

@@ -1504,6 +1504,16 @@ TEST(TriangleMesh, IsSelfIntersecting) {
                        {0.1, 0.1, 0}, {0.1, 1.1, 0}, {1.1, 0.1, 0}};
     mesh1.triangles_ = {{0, 1, 2}, {3, 4, 5}};
     EXPECT_TRUE(mesh1.IsSelfIntersecting());
+
+    // Regression for #5117: the AABBs touch at x = 1, but the triangles do
+    // not. The narrow phase must reject this pair.
+    geometry::TriangleMesh mesh2;
+    mesh2.vertices_ = {{0.0, 0.13918686, 1.0}, {0.0, 0.0, 1.1270161},
+                       {1.0, 0.0, 1.0284119},  {1.0, 1.1269569, 0.0},
+                       {1.0, 0.03113556, 1.0}, {2.0, 1.0189056, 0.0}};
+    mesh2.triangles_ = {{0, 1, 2}, {3, 4, 5}};
+    EXPECT_FALSE(mesh2.IsSelfIntersecting());
+    EXPECT_TRUE(mesh2.GetSelfIntersectingTriangles().empty());
 }
 
 TEST(TriangleMesh, GetVolume) {

python/test/geometry/test_trianglemesh.py

@@ -0,0 +1,32 @@
+# ----------------------------------------------------------------------------
+# -                        Open3D: www.open3d.org                            -
+# ----------------------------------------------------------------------------
+# Copyright (c) 2018-2024 www.open3d.org
+# SPDX-License-Identifier: MIT
+# ----------------------------------------------------------------------------
+
+import math
+
+import numpy as np
+import open3d as o3d
+
+
+def test_self_intersection_issue_5117():
+    vertices = np.array([[0.0, 0.13918686, 1.0], [0.0, 0.0, 1.1270161],
+                         [1.0, 0.0, 1.0284119], [1.0, 1.1269569, 0.0],
+                         [1.0, 0.03113556, 1.0], [2.0, 1.0189056, 0.0]])
+    triangles = np.array([[0, 1, 2], [3, 4, 5]])
+    mesh = o3d.geometry.TriangleMesh(o3d.utility.Vector3dVector(vertices),
+                                     o3d.utility.Vector3iVector(triangles))
+
+    # Classification of the separated issue #5117 triangles must remain false
+    # before and after a rigid rotation.
+    assert not mesh.is_self_intersecting()
+    assert len(mesh.get_self_intersecting_triangles()) == 0
+
+    sqrt_half = math.sqrt(0.5)
+    rotation = np.array([[1.0, 0.0, 0.0], [0.0, sqrt_half, sqrt_half],
+                         [0.0, -sqrt_half, sqrt_half]])
+    mesh.rotate(rotation)
+    assert not mesh.is_self_intersecting()
+    assert len(mesh.get_self_intersecting_triangles()) == 0