Fix T82301, exact boolean fail on cube with bump.

The code for determining coplanar clusters had a bug where it would
miss some triangles. The fix for now is to just put triangles in
the cluster if their bounding boxes overlap. This works but maybe
makes clusters bigger then they have to be. I'll follow this commit
with work on making the CDT routine faster when using exact arithmetic.
Also removed a lot of unused code, and added some new intersect
performance tests.

1 files changed, 2 insertions, 442 deletions

diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index a777833dff4..5f7258ebb6a 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -1107,254 +1107,6 @@ static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis)
 }
 
 /**
-   Is a point in the interior of a 2d triangle or on one of its
- * edges but not either endpoint of the edge?
- * orient[pi][i] is the orientation test of the point pi against
- * the side of the triangle starting at index i.
- * Assume the triangle is non-degenerate and CCW-oriented.
- * Then answer is true if p is left of or on all three of triangle a's edges,
- * and strictly left of at least on of them.
- */
-static bool non_trivially_2d_point_in_tri(const int orients[3][3], int pi)
-{
-  int p_left_01 = orients[pi][0];
-  int p_left_12 = orients[pi][1];
-  int p_left_20 = orients[pi][2];
-  return (p_left_01 >= 0 && p_left_12 >= 0 && p_left_20 >= 0 &&
-          (p_left_01 + p_left_12 + p_left_20) >= 2);
-}
-
-/**
- * Given orients as defined in non_trivially_2d_intersect, do the triangles
- * overlap in a "hex" pattern? That is, the overlap region is a hexagon, which
- * one gets by having, each point of one triangle being strictly right-of one
- * edge of the other and strictly left of the other two edges; and vice versa.
- * In addition, it must not be the case that all of the points of one triangle
- * are totally on the outside of one edge of the other triangle, and vice versa.
- */
-static bool non_trivially_2d_hex_overlap(int orients[2][3][3])
-{
-  for (int ab = 0; ab < 2; ++ab) {
-    for (int i = 0; i < 3; ++i) {
-      bool ok = orients[ab][i][0] + orients[ab][i][1] + orients[ab][i][2] == 1 &&
-                orients[ab][i][0] != 0 && orients[ab][i][1] != 0 && orients[i][2] != 0;
-      if (!ok) {
-        return false;
-      }
-      int s = orients[ab][0][i] + orients[ab][1][i] + orients[ab][2][i];
-      if (s == -3) {
-        return false;
-      }
-    }
-  }
-  return true;
-}
-
-/**
- * Given orients as defined in non_trivially_2d_intersect, do the triangles
- * have one shared edge in a "folded-over" configuration?
- * As well as a shared edge, the third vertex of one triangle needs to be
- * right-of one and left-of the other two edges of the other triangle.
- */
-static bool non_trivially_2d_shared_edge_overlap(int orients[2][3][3],
-                                                 const mpq2 *a[3],
-                                                 const mpq2 *b[3])
-{
-  for (int i = 0; i < 3; ++i) {
-    int in = (i + 1) % 3;
-    int inn = (i + 2) % 3;
-    for (int j = 0; j < 3; ++j) {
-      int jn = (j + 1) % 3;
-      int jnn = (j + 2) % 3;
-      if (*a[i] == *b[j] && *a[in] == *b[jn]) {
-        /* Edge from a[i] is shared with edge from b[j]. */
-        /* See if a[inn] is right-of or on one of the other edges of b.
-         * If it is on, then it has to be right-of or left-of the shared edge,
-         * depending on which edge it is. */
-        if (orients[0][inn][jn] < 0 || orients[0][inn][jnn] < 0) {
-          return true;
-        }
-        if (orients[0][inn][jn] == 0 && orients[0][inn][j] == 1) {
-          return true;
-        }
-        if (orients[0][inn][jnn] == 0 && orients[0][inn][j] == -1) {
-          return true;
-        }
-        /* Similarly for `b[jnn]`. */
-        if (orients[1][jnn][in] < 0 || orients[1][jnn][inn] < 0) {
-          return true;
-        }
-        if (orients[1][jnn][in] == 0 && orients[1][jnn][i] == 1) {
-          return true;
-        }
-        if (orients[1][jnn][inn] == 0 && orients[1][jnn][i] == -1) {
-          return true;
-        }
-      }
-    }
-  }
-  return false;
-}
-
-/**
- * Are the triangles the same, perhaps with some permutation of vertices?
- */
-static bool same_triangles(const mpq2 *a[3], const mpq2 *b[3])
-{
-  for (int i = 0; i < 3; ++i) {
-    if (a[0] == b[i] && a[1] == b[(i + 1) % 3] && a[2] == b[(i + 2) % 3]) {
-      return true;
-    }
-  }
-  return false;
-}
-
-/**
- * Do 2d triangles (a[0], a[1], a[2]) and (b[0], b[1], b2[2]) intersect at more than just shared
- * vertices or a shared edge? This is true if any point of one triangle is non-trivially inside the
- * other. NO: that isn't quite sufficient: there is also the case where the verts are all mutually
- * outside the other's triangle, but there is a hexagonal overlap region where they overlap.
- */
-static bool non_trivially_2d_intersect(const mpq2 *a[3], const mpq2 *b[3])
-{
-  /* TODO: Could experiment with trying bounding box tests before these.
-   * TODO: Find a less expensive way than 18 orient tests to do this. */
-
-  /* `orients[0][ai][bi]` is orient of point `a[ai]` compared to segment starting at `b[bi]`.
-   * `orients[1][bi][ai]` is orient of point `b[bi]` compared to segment starting at `a[ai]`. */
-  int orients[2][3][3];
-  for (int ab = 0; ab < 2; ++ab) {
-    for (int ai = 0; ai < 3; ++ai) {
-      for (int bi = 0; bi < 3; ++bi) {
-        if (ab == 0) {
-          orients[0][ai][bi] = orient2d(*b[bi], *b[(bi + 1) % 3], *a[ai]);
-        }
-        else {
-          orients[1][bi][ai] = orient2d(*a[ai], *a[(ai + 1) % 3], *b[bi]);
-        }
-      }
-    }
-  }
-  return non_trivially_2d_point_in_tri(orients[0], 0) ||
-         non_trivially_2d_point_in_tri(orients[0], 1) ||
-         non_trivially_2d_point_in_tri(orients[0], 2) ||
-         non_trivially_2d_point_in_tri(orients[1], 0) ||
-         non_trivially_2d_point_in_tri(orients[1], 1) ||
-         non_trivially_2d_point_in_tri(orients[1], 2) || non_trivially_2d_hex_overlap(orients) ||
-         non_trivially_2d_shared_edge_overlap(orients, a, b) || same_triangles(a, b);
-  return true;
-}
-
-/**
- * Does triangle t in tm non-trivially non-co-planar intersect any triangle
- * in `CoplanarCluster cl`? Assume t is known to be in the same plane as all
- * the triangles in cl, and that proj_axis is a good axis to project down
- * to solve this problem in 2d.
- */
-static bool non_trivially_coplanar_intersects(const IMesh &tm,
-                                              int t,
-                                              const CoplanarCluster &cl,
-                                              int proj_axis,
-                                              const Map<std::pair<int, int>, ITT_value> &itt_map)
-{
-  const Face &tri = *tm.face(t);
-  mpq2 v0 = project_3d_to_2d(tri[0]->co_exact, proj_axis);
-  mpq2 v1 = project_3d_to_2d(tri[1]->co_exact, proj_axis);
-  mpq2 v2 = project_3d_to_2d(tri[2]->co_exact, proj_axis);
-  if (orient2d(v0, v1, v2) != 1) {
-    mpq2 tmp = v1;
-    v1 = v2;
-    v2 = tmp;
-  }
-  for (const int cl_t : cl) {
-    if (!itt_map.contains(std::pair<int, int>(t, cl_t)) &&
-        !itt_map.contains(std::pair<int, int>(cl_t, t))) {
-      continue;
-    }
-    const Face &cl_tri = *tm.face(cl_t);
-    mpq2 ctv0 = project_3d_to_2d(cl_tri[0]->co_exact, proj_axis);
-    mpq2 ctv1 = project_3d_to_2d(cl_tri[1]->co_exact, proj_axis);
-    mpq2 ctv2 = project_3d_to_2d(cl_tri[2]->co_exact, proj_axis);
-    if (orient2d(ctv0, ctv1, ctv2) != 1) {
-      mpq2 tmp = ctv1;
-      ctv1 = ctv2;
-      ctv2 = tmp;
-    }
-    const mpq2 *v[] = {&v0, &v1, &v2};
-    const mpq2 *ctv[] = {&ctv0, &ctv1, &ctv2};
-    if (non_trivially_2d_intersect(v, ctv)) {
-      return true;
-    }
-  }
-  return false;
-}
-
-/* Keeping this code for a while, but for now, almost all
- * trivial intersects are found before calling intersect_tri_tri now.
- */
-#  if 0
-/**
- * Do tri1 and tri2 intersect at all, and if so, is the intersection
- * something other than a common vertex or a common edge?
- * The \a itt value is the result of calling intersect_tri_tri on tri1, tri2.
- */
-static bool non_trivial_intersect(const ITT_value &itt, const Face * tri1, const Face * tri2)
-{
-  if (itt.kind == INONE) {
-    return false;
-  }
-  const Face * tris[2] = {tri1, tri2};
-  if (itt.kind == IPOINT) {
-    bool has_p_as_vert[2] {false, false};
-    for (int i = 0; i < 2; ++i) {
-      for (const Vert * v : *tris[i]) {
-        if (itt.p1 == v->co_exact) {
-          has_p_as_vert[i] = true;
-          break;
-        }
-      }
-    }
-    return !(has_p_as_vert[0] && has_p_as_vert[1]);
-  }
-  if (itt.kind == ISEGMENT) {
-    bool has_seg_as_edge[2] = {false, false};
-    for (int i = 0; i < 2; ++i) {
-      const Face &t = *tris[i];
-      for (int pos : t.index_range()) {
-        int nextpos = t.next_pos(pos);
-        if ((itt.p1 == t[pos]->co_exact && itt.p2 == t[nextpos]->co_exact) ||
-            (itt.p2 == t[pos]->co_exact && itt.p1 == t[nextpos]->co_exact)) {
-          has_seg_as_edge[i] = true;
-          break;
-        }
-      }
-    }
-    return !(has_seg_as_edge[0] && has_seg_as_edge[1]);
-  }
-  BLI_assert(itt.kind == ICOPLANAR);
-  /* TODO: refactor this common code with code above. */
-  int proj_axis = mpq3::dominant_axis(tri1->plane.norm_exact);
-  mpq2 tri_2d[2][3];
-  for (int i = 0; i < 2; ++i) {
-    mpq2 v0 = project_3d_to_2d((*tris[i])[0]->co_exact, proj_axis);
-    mpq2 v1 = project_3d_to_2d((*tris[i])[1]->co_exact, proj_axis);
-    mpq2 v2 = project_3d_to_2d((*tris[i])[2]->co_exact, proj_axis);
-    if (mpq2::orient2d(v0, v1, v2) != 1) {
-      mpq2 tmp = v1;
-      v1 = v2;
-      v2 = tmp;
-    }
-    tri_2d[i][0] = v0;
-    tri_2d[i][1] = v1;
-    tri_2d[i][2] = v2;
-  }
-  const mpq2 *va[] = {&tri_2d[0][0], &tri_2d[0][1], &tri_2d[0][2]};
-  const mpq2 *vb[] = {&tri_2d[1][0], &tri_2d[1][1], &tri_2d[1][2]};
-  return non_trivially_2d_intersect(va, vb);
-}
-#  endif
-
-/**
  * The sup and index functions are defined in the paper:
  * EXACT GEOMETRIC COMPUTATION USING CASCADING, by
  * Burnikel, Funke, and Seel. They are used to find absolute
@@ -1414,119 +1166,6 @@ constexpr int index_dot_plane_coords = 15;
  */
 constexpr int index_dot_cross = 11;
 
-/* Not using this at the moment. Leaving it for a bit in case we want it again. */
-#  if 0
-static double supremum_dot(const double3 &a, const double3 &b)
-{
-  double3 abs_a = double3::abs(a);
-  double3 abs_b = double3::abs(b);
-  return double3::dot(abs_a, abs_b);
-}
-
-static double supremum_orient3d(const double3 &a,
-                                const double3 &b,
-                                const double3 &c,
-                                const double3 &d)
-{
-  double3 abs_a = double3::abs(a);
-  double3 abs_b = double3::abs(b);
-  double3 abs_c = double3::abs(c);
-  double3 abs_d = double3::abs(d);
-  double adx = abs_a[0] + abs_d[0];
-  double bdx = abs_b[0] + abs_d[0];
-  double cdx = abs_c[0] + abs_d[0];
-  double ady = abs_a[1] + abs_d[1];
-  double bdy = abs_b[1] + abs_d[1];
-  double cdy = abs_c[1] + abs_d[1];
-  double adz = abs_a[2] + abs_d[2];
-  double bdz = abs_b[2] + abs_d[2];
-  double cdz = abs_c[2] + abs_d[2];
-
-  double bdxcdy = bdx * cdy;
-  double cdxbdy = cdx * bdy;
-
-  double cdxady = cdx * ady;
-  double adxcdy = adx * cdy;
-
-  double adxbdy = adx * bdy;
-  double bdxady = bdx * ady;
-
-  double det = adz * (bdxcdy + cdxbdy) + bdz * (cdxady + adxcdy) + cdz * (adxbdy + bdxady);
-  return det;
-}
-
-/** Actually index_orient3d = 10 + 4 * (max degree of input coordinates) */
-constexpr int index_orient3d = 14;
-
-/**
- * Return the approximate orient3d of the four double3's, with
- * the guarantee that if the value is -1 or 1 then the underlying
- * mpq3 test would also have returned that value.
- * When the return value is 0, we are not sure of the sign.
- */
-static int filter_orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
-{
-  double o3dfast = orient3d_fast(a, b, c, d);
-  if (o3dfast == 0.0) {
-    return 0;
-  }
-  double err_bound = supremum_orient3d(a, b, c, d) * index_orient3d * DBL_EPSILON;
-  if (fabs(o3dfast) > err_bound) {
-    return o3dfast > 0.0 ? 1 : -1;
-  }
-  return 0;
-}
-
-/**
- * Return the approximate orient3d of the triangle plane points and v, with
- * the guarantee that if the value is -1 or 1 then the underlying
- * mpq3 test would also have returned that value.
- * When the return value is 0, we are not sure of the sign.
- */
-static int filter_tri_plane_vert_orient3d(const Face &tri, const Vert *v)
-{
-  return filter_orient3d(tri[0]->co, tri[1]->co, tri[2]->co, v->co);
-}
-
-/**
- * Are vectors a and b parallel or nearly parallel?
- * This routine should only return false if we are certain
- * that they are not parallel, taking into account the
- * possible numeric errors and input value approximation.
- */
-static bool near_parallel_vecs(const double3 &a, const double3 &b)
-{
-  double3 cr = double3::cross_high_precision(a, b);
-  double cr_len_sq = cr.length_squared();
-  if (cr_len_sq == 0.0) {
-    return true;
-  }
-  double err_bound = supremum_dot_cross(a, b) * index_dot_cross * DBL_EPSILON;
-  if (cr_len_sq > err_bound) {
-    return false;
-  }
-  return true;
-}
-
-/**
- * Return true if we are sure that dot(a,b) > 0, taking into
- * account the error bounds due to numeric errors and input value
- * approximation.
- */
-static bool dot_must_be_positive(const double3 &a, const double3 &b)
-{
-  double d = double3::dot(a, b);
-  if (d <= 0.0) {
-    return false;
-  }
-  double err_bound = supremum_dot(a, b) * index_dot_plane_coords * DBL_EPSILON;
-  if (d > err_bound) {
-    return true;
-  }
-  return false;
-}
-#  endif
-
 /**
  * Return the approximate side of point p on a plane with normal plane_no and point plane_p.
  * The answer will be 1 if p is definitely above the plane, -1 if it is definitely below.
@@ -1556,83 +1195,6 @@ static int filter_plane_side(const double3 &p,
   return 0;
 }
 
-/* Not using this at the moment. Leave it here for a while in case we want it again. */
-#  if 0
-/**
- * A fast, non-exhaustive test for non_trivial intersection.
- * If this returns false then we are sure that tri1 and tri2
- * do not intersect. If it returns true, they may or may not
- * non-trivially intersect.
- * We assume that bounding box overlap tests have already been
- * done, so don't repeat those here. This routine is checking
- * for the very common cases (when doing mesh self-intersect)
- * where triangles share an edge or a vertex, but don't
- * otherwise intersect.
- */
-static bool may_non_trivially_intersect(Face *t1, Face *t2)
-{
-  Face &tri1 = *t1;
-  Face &tri2 = *t2;
-  BLI_assert(t1->plane_populated() && t2->plane_populated());
-  Face::FacePos share1_pos[3];
-  Face::FacePos share2_pos[3];
-  int n_shared = 0;
-  for (Face::FacePos p1 = 0; p1 < 3; ++p1) {
-    const Vert *v1 = tri1[p1];
-    for (Face::FacePos p2 = 0; p2 < 3; ++p2) {
-      const Vert *v2 = tri2[p2];
-      if (v1 == v2) {
-        share1_pos[n_shared] = p1;
-        share2_pos[n_shared] = p2;
-        ++n_shared;
-      }
-    }
-  }
-  if (n_shared == 2) {
-    /* t1 and t2 share an entire edge.
-     * If their normals are not parallel, they cannot non-trivially intersect. */
-    if (!near_parallel_vecs(tri1.plane->norm, tri2.plane->norm)) {
-      return false;
-    }
-    /* The normals are parallel or nearly parallel.
-     * If the normals are in the same direction and the edges have opposite
-     * directions in the two triangles, they cannot non-trivially intersect. */
-    bool erev1 = tri1.prev_pos(share1_pos[0]) == share1_pos[1];
-    bool erev2 = tri2.prev_pos(share2_pos[0]) == share2_pos[1];
-    if (erev1 != erev2) {
-      if (dot_must_be_positive(tri1.plane->norm, tri2.plane->norm)) {
-        return false;
-      }
-    }
-  }
-  else if (n_shared == 1) {
-    /* t1 and t2 share a vertex, but not an entire edge.
-     * If the two non-shared verts of t2 are both on the same
-     * side of tri1's plane, then they cannot non-trivially intersect.
-     * (There are some other cases that could be caught here but
-     * they are more expensive to check). */
-    Face::FacePos p = share2_pos[0];
-    const Vert *v2a = p == 0 ? tri2[1] : tri2[0];
-    const Vert *v2b = (p == 0 || p == 1) ? tri2[2] : tri2[1];
-    int o1 = filter_tri_plane_vert_orient3d(tri1, v2a);
-    int o2 = filter_tri_plane_vert_orient3d(tri1, v2b);
-    if (o1 == o2 && o1 != 0) {
-      return false;
-    }
-    p = share1_pos[0];
-    const Vert *v1a = p == 0 ? tri1[1] : tri1[0];
-    const Vert *v1b = (p == 0 || p == 1) ? tri1[2] : tri1[1];
-    o1 = filter_tri_plane_vert_orient3d(tri2, v1a);
-    o2 = filter_tri_plane_vert_orient3d(tri2, v1b);
-    if (o1 == o2 && o1 != 0) {
-      return false;
-    }
-  }
-  /* We weren't able to prove that any intersection is trivial. */
-  return true;
-}
-#  endif
-
 /*
  * interesect_tri_tri and helper functions.
  * This code uses the algorithm of Guigue and Devillers, as described
@@ -2852,7 +2414,6 @@ static CoplanarClusterInfo find_clusters(const IMesh &tm,
       if (dbg_level > 0) {
         std::cout << "already has " << curcls.size() << " clusters\n";
       }
-      int proj_axis = mpq3::dominant_axis(tplane.norm_exact);
       /* Partition `curcls` into those that intersect t non-trivially, and those that don't. */
       Vector<CoplanarCluster *> int_cls;
       Vector<CoplanarCluster *> no_int_cls;
@@ -2860,8 +2421,7 @@ static CoplanarClusterInfo find_clusters(const IMesh &tm,
         if (dbg_level > 1) {
           std::cout << "consider intersecting with cluster " << cl << "\n";
         }
-        if (bbs_might_intersect(tri_bb[t], cl.bounding_box()) &&
-            non_trivially_coplanar_intersects(tm, t, cl, proj_axis, itt_map)) {
+        if (bbs_might_intersect(tri_bb[t], cl.bounding_box())) {
           if (dbg_level > 1) {
             std::cout << "append to int_cls\n";
           }
@@ -3317,6 +2877,6 @@ static void dump_perfdata(void)
 }
 #  endif
 
-};  // namespace blender::meshintersect
+}  // namespace blender::meshintersect
 
 #endif  // WITH_GMP