Fix T87554 Exact Boolean performance bug.

There was a quadratic algorithm extracting triangles from a coplanar
cluster. This is now linear.
Also found and fixed a bug in the same area related to the triangulator
added recently: it didn't get the right correspondence between new
edges and original edges.

1 files changed, 391 insertions, 92 deletions

diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index ce3a5b55f98..3558dfad81d 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -39,6 +39,7 @@
 #  include "BLI_math_mpq.hh"
 #  include "BLI_mpq2.hh"
 #  include "BLI_mpq3.hh"
+#  include "BLI_polyfill_2d.h"
 #  include "BLI_set.hh"
 #  include "BLI_span.hh"
 #  include "BLI_task.h"
@@ -1576,6 +1577,9 @@ struct CDT_data {
   Vector<bool> is_reversed;
   /** Result of running CDT on input with (vert, edge, face). */
   CDT_result<mpq_class> cdt_out;
+  /** To speed up get_cdt_edge_orig, sometimes populate this map from vertex pair to output edge.
+   */
+  Map<std::pair<int, int>, int> verts_to_edge;
   int proj_axis;
 };
 
@@ -1734,6 +1738,32 @@ static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm,
   return ans;
 }
 
+/* Return a copy of the argument with the integers ordered in ascending order. */
+static inline std::pair<int, int> sorted_int_pair(std::pair<int, int> pair)
+{
+  if (pair.first <= pair.second) {
+    return pair;
+  }
+  else {
+    return std::pair<int, int>(pair.second, pair.first);
+  }
+}
+
+/**
+ * Build cd.verts_to_edge to map from a pair of cdt output indices to an index in cd.cdt_out.edge.
+ * Order the vertex indices so that the smaller one is first in the pair.
+ */
+static void populate_cdt_edge_map(Map<std::pair<int, int>, int> &verts_to_edge,
+                                  const CDT_result<mpq_class> &cdt_out)
+{
+  verts_to_edge.reserve(cdt_out.edge.size());
+  for (int e : cdt_out.edge.index_range()) {
+    std::pair<int, int> vpair = sorted_int_pair(cdt_out.edge[e]);
+    /* There should be only one edge for each vertex pair. */
+    verts_to_edge.add(vpair, e);
+  }
+}
+
 /**
  * Fills in cd.cdt_out with result of doing the cdt calculation on (vert, edge, face).
  */
@@ -1766,6 +1796,10 @@ static void do_cdt(CDT_data &cd)
   }
   cdt_in.epsilon = 0; /* TODO: needs attention for non-exact T. */
   cd.cdt_out = blender::meshintersect::delaunay_2d_calc(cdt_in, CDT_INSIDE);
+  constexpr int make_edge_map_threshold = 15;
+  if (cd.cdt_out.edge.size() >= make_edge_map_threshold) {
+    populate_cdt_edge_map(cd.verts_to_edge, cd.cdt_out);
+  }
   if (dbg_level > 0) {
     std::cout << "\nCDT result\nVerts:\n";
     for (int i : cd.cdt_out.vert.index_range()) {
@@ -1802,39 +1836,52 @@ static int get_cdt_edge_orig(
 {
   int foff = cd.cdt_out.face_edge_offset;
   *r_is_intersect = false;
-  for (int e : cd.cdt_out.edge.index_range()) {
-    std::pair<int, int> edge = cd.cdt_out.edge[e];
-    if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) {
-      /* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */
-      /* TODO: if edge has origs from more than on part of the nary input,
-       * then want to set *r_is_intersect to true. */
-      for (int orig_index : cd.cdt_out.edge_orig[e]) {
-        /* orig_index encodes the triangle and pos within the triangle of the input edge. */
-        if (orig_index >= foff) {
-          int in_face_index = (orig_index / foff) - 1;
-          int pos = orig_index % foff;
-          /* We need to retrieve the edge orig field from the Face used to populate the
-           * in_face_index'th face of the CDT, at the pos'th position of the face. */
-          int in_tm_face_index = cd.input_face[in_face_index];
-          BLI_assert(in_tm_face_index < in_tm.face_size());
-          const Face *facep = in_tm.face(in_tm_face_index);
-          BLI_assert(pos < facep->size());
-          bool is_rev = cd.is_reversed[in_face_index];
-          int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos];
-          if (eorig != NO_INDEX) {
-            return eorig;
-          }
-        }
-        else {
-          /* This edge came from an edge input to the CDT problem,
-           * so it is an intersect edge. */
-          *r_is_intersect = true;
-          /* TODO: maybe there is an orig index:
-           * This happens if an input edge was formed by an input face having
-           * an edge that is co-planar with the cluster, while the face as a whole is not. */
-          return NO_INDEX;
-        }
+  int e = NO_INDEX;
+  if (cd.verts_to_edge.size() > 0) {
+    /* Use the populated map to find the edge, if any, between vertices i0 and i1. */
+    std::pair<int, int> vpair = sorted_int_pair(std::pair<int, int>(i0, i1));
+    e = cd.verts_to_edge.lookup_default(vpair, NO_INDEX);
+  }
+  else {
+    for (int ee : cd.cdt_out.edge.index_range()) {
+      std::pair<int, int> edge = cd.cdt_out.edge[ee];
+      if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) {
+        e = ee;
+        break;
       }
+    }
+  }
+  if (e == NO_INDEX) {
+    return NO_INDEX;
+  }
+
+  /* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */
+  /* TODO: if edge has origs from more than on part of the nary input,
+   * then want to set *r_is_intersect to true. */
+  for (int orig_index : cd.cdt_out.edge_orig[e]) {
+    /* orig_index encodes the triangle and pos within the triangle of the input edge. */
+    if (orig_index >= foff) {
+      int in_face_index = (orig_index / foff) - 1;
+      int pos = orig_index % foff;
+      /* We need to retrieve the edge orig field from the Face used to populate the
+       * in_face_index'th face of the CDT, at the pos'th position of the face. */
+      int in_tm_face_index = cd.input_face[in_face_index];
+      BLI_assert(in_tm_face_index < in_tm.face_size());
+      const Face *facep = in_tm.face(in_tm_face_index);
+      BLI_assert(pos < facep->size());
+      bool is_rev = cd.is_reversed[in_face_index];
+      int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos];
+      if (eorig != NO_INDEX) {
+        return eorig;
+      }
+    }
+    else {
+      /* This edge came from an edge input to the CDT problem,
+       * so it is an intersect edge. */
+      *r_is_intersect = true;
+      /* TODO: maybe there is an orig index:
+       * This happens if an input edge was formed by an input face having
+       * an edge that is co-planar with the cluster, while the face as a whole is not. */
       return NO_INDEX;
     }
   }
@@ -1842,6 +1889,256 @@ static int get_cdt_edge_orig(
 }
 
 /**
+ * Make a Face from the CDT output triangle cdt_out_t, which has corresponding input triangle
+ * cdt_in_t. The triangle indices that are in cd.input_face are from IMesh tm.
+ * Return a pointer to the newly constructed Face.
+ */
+static Face *cdt_tri_as_imesh_face(
+    int cdt_out_t, int cdt_in_t, const CDT_data &cd, const IMesh &tm, IMeshArena *arena)
+{
+  const CDT_result<mpq_class> &cdt_out = cd.cdt_out;
+  int t_orig = tm.face(cd.input_face[cdt_in_t])->orig;
+  BLI_assert(cdt_out.face[cdt_out_t].size() == 3);
+  int i0 = cdt_out.face[cdt_out_t][0];
+  int i1 = cdt_out.face[cdt_out_t][1];
+  int i2 = cdt_out.face[cdt_out_t][2];
+  mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]);
+  mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]);
+  mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]);
+  /* No need to provide an original index: if coord matches
+   * an original one, then it will already be in the arena
+   * with the correct orig field. */
+  const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX);
+  const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX);
+  const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX);
+  Face *facep;
+  bool is_isect0;
+  bool is_isect1;
+  bool is_isect2;
+  if (cd.is_reversed[cdt_in_t]) {
+    int oe0 = get_cdt_edge_orig(i0, i2, cd, tm, &is_isect0);
+    int oe1 = get_cdt_edge_orig(i2, i1, cd, tm, &is_isect1);
+    int oe2 = get_cdt_edge_orig(i1, i0, cd, tm, &is_isect2);
+    facep = arena->add_face(
+        {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
+  }
+  else {
+    int oe0 = get_cdt_edge_orig(i0, i1, cd, tm, &is_isect0);
+    int oe1 = get_cdt_edge_orig(i1, i2, cd, tm, &is_isect1);
+    int oe2 = get_cdt_edge_orig(i2, i0, cd, tm, &is_isect2);
+    facep = arena->add_face(
+        {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
+  }
+  facep->populate_plane(false);
+  return facep;
+}
+
+/**
+ * Tessellate face f into triangles and return an array of `const Face *`
+ * giving that triangulation. Intended to be used when f has > 4 vertices.
+ * Care is taken so that the original edge index associated with
+ * each edge in the output triangles either matches the original edge
+ * for the (identical) edge of f, or else is -1. So diagonals added
+ * for triangulation can later be identified by having #NO_INDEX for original.
+ *
+ * This code uses Blenlib's BLI_polyfill_calc to do triangulation, and is therefore quite fast.
+ * Unfortunately, it can product degenerate triangles that mesh_intersect will remove, leaving
+ * the mesh non-PWN.
+ */
+static Array<Face *> polyfill_triangulate_poly(Face *f, IMeshArena *arena)
+{
+  /* Similar to loop body in BM_mesh_calc_tesselation. */
+  int flen = f->size();
+  BLI_assert(flen > 4);
+  if (!f->plane_populated()) {
+    f->populate_plane(false);
+  }
+  /* Project along negative face normal so (x,y) can be used in 2d. */
+  const double3 &poly_normal = f->plane->norm;
+  float no[3] = {float(poly_normal[0]), float(poly_normal[1]), float(poly_normal[2])};
+  normalize_v3(no);
+  float axis_mat[3][3];
+  float(*projverts)[2];
+  unsigned int(*tris)[3];
+  const int totfilltri = flen - 2;
+  /* Prepare projected vertices and array to receive triangles in tesselation. */
+  tris = static_cast<unsigned int(*)[3]>(MEM_malloc_arrayN(totfilltri, sizeof(*tris), __func__));
+  projverts = static_cast<float(*)[2]>(MEM_malloc_arrayN(flen, sizeof(*projverts), __func__));
+  axis_dominant_v3_to_m3_negate(axis_mat, no);
+  for (int j = 0; j < flen; ++j) {
+    const double3 &dco = (*f)[j]->co;
+    float co[3] = {float(dco[0]), float(dco[1]), float(dco[2])};
+    mul_v2_m3v3(projverts[j], axis_mat, co);
+  }
+  BLI_polyfill_calc(projverts, flen, 1, tris);
+  /* Put tesselation triangles into Face form. Record original edges where they exist. */
+  Array<Face *> ans(totfilltri);
+  for (int t = 0; t < totfilltri; ++t) {
+    unsigned int *tri = tris[t];
+    int eo[3];
+    const Vert *v[3];
+    for (int k = 0; k < 3; k++) {
+      BLI_assert(tri[k] < flen);
+      v[k] = (*f)[tri[k]];
+      /* If tri edge goes between two successive indices in
+       * the original face, then it is an original edge. */
+      if ((tri[k] + 1) % flen == tri[(k + 1) % 3]) {
+        eo[k] = f->edge_orig[tri[k]];
+      }
+      else {
+        eo[k] = NO_INDEX;
+      }
+      ans[t] = arena->add_face(
+          {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false});
+    }
+  }
+
+  MEM_freeN(tris);
+  MEM_freeN(projverts);
+
+  return ans;
+}
+
+/**
+ * Tessellate face f into triangles and return an array of `const Face *`
+ * giving that triangulation.
+ * Care is taken so that the original edge index associated with
+ * each edge in the output triangles either matches the original edge
+ * for the (identical) edge of f, or else is -1. So diagonals added
+ * for triangulation can later be identified by having #NO_INDEX for original.
+ *
+ * The method used is to use the CDT triangulation. Usually that triangulation
+ * will only use the existing vertices. However, if the face self-intersects
+ * then the CDT triangulation will include the intersection points.
+ * If this happens, we use the polyfill triangulator instead. We don't
+ * use the polyfill triangulator by default because it can create degenerate
+ * triangles (which we can handle but they'll create non-manifold meshes).
+ */
+static Array<Face *> triangulate_poly(Face *f, IMeshArena *arena)
+{
+  int flen = f->size();
+  CDT_input<mpq_class> cdt_in;
+  cdt_in.vert = Array<mpq2>(flen);
+  cdt_in.face = Array<Vector<int>>(1);
+  cdt_in.face[0].reserve(flen);
+  for (int i : f->index_range()) {
+    cdt_in.face[0].append(i);
+  }
+  /* Project poly along dominant axis of normal to get 2d coords. */
+  if (!f->plane_populated()) {
+    f->populate_plane(false);
+  }
+  const double3 &poly_normal = f->plane->norm;
+  int axis = double3::dominant_axis(poly_normal);
+  /* If project down y axis as opposed to x or z, the orientation
+   * of the polygon will be reversed.
+   * Yet another reversal happens if the poly normal in the dominant
+   * direction is opposite that of the positive dominant axis. */
+  bool rev1 = (axis == 1);
+  bool rev2 = poly_normal[axis] < 0;
+  bool rev = rev1 ^ rev2;
+  for (int i = 0; i < flen; ++i) {
+    int ii = rev ? flen - i - 1 : i;
+    mpq2 &p2d = cdt_in.vert[ii];
+    int k = 0;
+    for (int j = 0; j < 3; ++j) {
+      if (j != axis) {
+        p2d[k++] = (*f)[ii]->co_exact[j];
+      }
+    }
+  }
+  CDT_result<mpq_class> cdt_out = delaunay_2d_calc(cdt_in, CDT_INSIDE);
+  int n_tris = cdt_out.face.size();
+  Array<Face *> ans(n_tris);
+  for (int t = 0; t < n_tris; ++t) {
+    int i_v_out[3];
+    const Vert *v[3];
+    int eo[3];
+    bool needs_steiner = false;
+    for (int i = 0; i < 3; ++i) {
+      i_v_out[i] = cdt_out.face[t][i];
+      if (cdt_out.vert_orig[i_v_out[i]].size() == 0) {
+        needs_steiner = true;
+        break;
+      }
+      v[i] = (*f)[cdt_out.vert_orig[i_v_out[i]][0]];
+    }
+    if (needs_steiner) {
+      /* Fall back on the polyfill triangulator. */
+      return polyfill_triangulate_poly(f, arena);
+    }
+    Map<std::pair<int, int>, int> verts_to_edge;
+    populate_cdt_edge_map(verts_to_edge, cdt_out);
+    int foff = cdt_out.face_edge_offset;
+    for (int i = 0; i < 3; ++i) {
+      std::pair<int, int> vpair(i_v_out[i], i_v_out[(i + 1) % 3]);
+      std::pair<int, int> vpair_canon = sorted_int_pair(vpair);
+      int e_out = verts_to_edge.lookup_default(vpair_canon, NO_INDEX);
+      BLI_assert(e_out != NO_INDEX);
+      eo[i] = NO_INDEX;
+      for (int orig : cdt_out.edge_orig[e_out]) {
+        if (orig >= foff) {
+          int pos = orig % foff;
+          BLI_assert(pos < f->size());
+          eo[i] = f->edge_orig[pos];
+          break;
+        }
+      }
+    }
+    if (rev) {
+      ans[t] = arena->add_face(
+          {v[0], v[2], v[1]}, f->orig, {eo[2], eo[1], eo[0]}, {false, false, false});
+    }
+    else {
+      ans[t] = arena->add_face(
+          {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false});
+    }
+  }
+  return ans;
+}
+
+/**
+ * Return an #IMesh that is a triangulation of a mesh with general
+ * polygonal faces, #IMesh.
+ * Added diagonals will be distinguishable by having edge original
+ * indices of #NO_INDEX.
+ */
+IMesh triangulate_polymesh(IMesh &imesh, IMeshArena *arena)
+{
+  Vector<Face *> face_tris;
+  constexpr int estimated_tris_per_face = 3;
+  face_tris.reserve(estimated_tris_per_face * imesh.face_size());
+  for (Face *f : imesh.faces()) {
+    /* Tessellate face f, following plan similar to #BM_face_calc_tesselation. */
+    int flen = f->size();
+    if (flen == 3) {
+      face_tris.append(f);
+    }
+    else if (flen == 4) {
+      const Vert *v0 = (*f)[0];
+      const Vert *v1 = (*f)[1];
+      const Vert *v2 = (*f)[2];
+      const Vert *v3 = (*f)[3];
+      int eo_01 = f->edge_orig[0];
+      int eo_12 = f->edge_orig[1];
+      int eo_23 = f->edge_orig[2];
+      int eo_30 = f->edge_orig[3];
+      Face *f0 = arena->add_face({v0, v1, v2}, f->orig, {eo_01, eo_12, -1}, {false, false, false});
+      Face *f1 = arena->add_face({v0, v2, v3}, f->orig, {-1, eo_23, eo_30}, {false, false, false});
+      face_tris.append(f0);
+      face_tris.append(f1);
+    }
+    else {
+      Array<Face *> tris = triangulate_poly(f, arena);
+      for (Face *tri : tris) {
+        face_tris.append(tri);
+      }
+    }
+  }
+  return IMesh(face_tris);
+}
+
+/**
  * Using the result of CDT in cd.cdt_out, extract an #IMesh representing the subdivision
  * of input triangle t, which should be an element of cd.input_face.
  */
@@ -1862,55 +2159,17 @@ static IMesh extract_subdivided_tri(const CDT_data &cd,
     BLI_assert(false);
     return IMesh();
   }
-  int t_orig = in_tm.face(t)->orig;
   constexpr int inline_buf_size = 20;
   Vector<Face *, inline_buf_size> faces;
   for (int f : cdt_out.face.index_range()) {
     if (cdt_out.face_orig[f].contains(t_in_cdt)) {
-      BLI_assert(cdt_out.face[f].size() == 3);
-      int i0 = cdt_out.face[f][0];
-      int i1 = cdt_out.face[f][1];
-      int i2 = cdt_out.face[f][2];
-      mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]);
-      mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]);
-      mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]);
-      /* No need to provide an original index: if coord matches
-       * an original one, then it will already be in the arena
-       * with the correct orig field. */
-      const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX);
-      const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX);
-      const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX);
-      Face *facep;
-      bool is_isect0;
-      bool is_isect1;
-      bool is_isect2;
-      if (cd.is_reversed[t_in_cdt]) {
-        int oe0 = get_cdt_edge_orig(i0, i2, cd, in_tm, &is_isect0);
-        int oe1 = get_cdt_edge_orig(i2, i1, cd, in_tm, &is_isect1);
-        int oe2 = get_cdt_edge_orig(i1, i0, cd, in_tm, &is_isect2);
-        facep = arena->add_face(
-            {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
-      }
-      else {
-        int oe0 = get_cdt_edge_orig(i0, i1, cd, in_tm, &is_isect0);
-        int oe1 = get_cdt_edge_orig(i1, i2, cd, in_tm, &is_isect1);
-        int oe2 = get_cdt_edge_orig(i2, i0, cd, in_tm, &is_isect2);
-        facep = arena->add_face(
-            {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2});
-      }
-      facep->populate_plane(false);
+      Face *facep = cdt_tri_as_imesh_face(f, t_in_cdt, cd, in_tm, arena);
       faces.append(facep);
     }
   }
   return IMesh(faces);
 }
 
-static IMesh extract_single_tri(const IMesh &tm, int t)
-{
-  Face *f = tm.face(t);
-  return IMesh({f});
-}
-
 static bool bvhtreeverlap_cmp(const BVHTreeOverlap &a, const BVHTreeOverlap &b)
 {
   if (a.indexA < b.indexA) {
@@ -2126,7 +2385,7 @@ static void calc_overlap_itts(Map<std::pair<int, int>, ITT_value> &itt_map,
 }
 
 /**
- * Data needed for parallelization of calc_subdivided_tris.
+ * Data needed for parallelization of calc_subdivided_non_cluster_tris.
  */
 struct OverlapTriRange {
   int tri_index;
@@ -2203,14 +2462,13 @@ static void calc_subdivided_tri_range_func(void *__restrict userdata,
  * r_tri_subdivided with the result of intersecting it with
  * all the other triangles in the mesh, if it intersects any others.
  * But don't do this for triangles that are part of a cluster.
- * Also, do nothing here if the answer is just the triangle itself.
  */
-static void calc_subdivided_tris(Array<IMesh> &r_tri_subdivided,
-                                 const IMesh &tm,
-                                 const Map<std::pair<int, int>, ITT_value> &itt_map,
-                                 const CoplanarClusterInfo &clinfo,
-                                 const TriOverlaps &ov,
-                                 IMeshArena *arena)
+static void calc_subdivided_non_cluster_tris(Array<IMesh> &r_tri_subdivided,
+                                             const IMesh &tm,
+                                             const Map<std::pair<int, int>, ITT_value> &itt_map,
+                                             const CoplanarClusterInfo &clinfo,
+                                             const TriOverlaps &ov,
+                                             IMeshArena *arena)
 {
   const int dbg_level = 0;
   if (dbg_level > 0) {
@@ -2250,6 +2508,54 @@ static void calc_subdivided_tris(Array<IMesh> &r_tri_subdivided,
   settings.use_threading = intersect_use_threading;
   BLI_task_parallel_range(
       0, overlap_tri_range_tot, &data, calc_subdivided_tri_range_func, &settings);
+  /* Now have to put in the triangles that are the same as the input ones, and not in clusters.
+   */
+  for (int t : tm.face_index_range()) {
+    if (r_tri_subdivided[t].face_size() == 0 && clinfo.tri_cluster(t) == NO_INDEX) {
+      r_tri_subdivided[t] = IMesh({tm.face(t)});
+    }
+  }
+}
+
+/**
+ * For each cluster in clinfo, extract the triangles from the cluster
+ * that correspond to each original triangle t that is part of the cluster,
+ * and put the resulting triangles into an IMesh in tri_subdivided[t].
+ * We have already done the CDT for the triangles in the cluster, whose
+ * result is in cluster_subdivided[c] for each cluster c.
+ */
+static void calc_cluster_tris(Array<IMesh> &tri_subdivided,
+                              const IMesh &tm,
+                              const CoplanarClusterInfo &clinfo,
+                              const Array<CDT_data> &cluster_subdivided,
+                              IMeshArena *arena)
+{
+  for (int c : clinfo.index_range()) {
+    const CoplanarCluster &cl = clinfo.cluster(c);
+    const CDT_data &cd = cluster_subdivided[c];
+    /* Each triangle in cluster c should be an input triangle in cd.input_faces.
+     * (See prepare_cdt_input_for_cluster.)
+     * So accumulate a Vector of Face* for each input face by going through the
+     * output faces and making a Face for each input face that it is part of.
+     * (The Boolean algorithm wants duplicates if a given output triangle is part
+     * of more than one input triangle.)
+     */
+    int n_cluster_tris = cl.tot_tri();
+    const CDT_result<mpq_class> &cdt_out = cd.cdt_out;
+    BLI_assert(cd.input_face.size() == n_cluster_tris);
+    Array<Vector<Face *>> face_vec(n_cluster_tris);
+    for (int cdt_out_t : cdt_out.face.index_range()) {
+      for (int cdt_in_t : cdt_out.face_orig[cdt_out_t]) {
+        Face *f = cdt_tri_as_imesh_face(cdt_out_t, cdt_in_t, cd, tm, arena);
+        face_vec[cdt_in_t].append(f);
+      }
+    }
+    for (int cdt_in_t : cd.input_face.index_range()) {
+      int tm_t = cd.input_face[cdt_in_t];
+      BLI_assert(tri_subdivided[tm_t].face_size() == 0);
+      tri_subdivided[tm_t] = IMesh(face_vec[cdt_in_t]);
+    }
+  }
 }
 
 static CDT_data calc_cluster_subdivided(const CoplanarClusterInfo &clinfo,
@@ -2622,10 +2928,11 @@ IMesh trimesh_nary_intersect(const IMesh &tm_in,
   doperfmax(2, tri_ov.overlap().size());
 #  endif
   Array<IMesh> tri_subdivided(tm_clean->face_size());
-  calc_subdivided_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena);
+  calc_subdivided_non_cluster_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena);
 #  ifdef PERFDEBUG
   double subdivided_tris_time = PIL_check_seconds_timer();
-  std::cout << "subdivided tris found, time = " << subdivided_tris_time - itt_time << "\n";
+  std::cout << "subdivided non-cluster tris found, time = " << subdivided_tris_time - itt_time
+            << "\n";
 #  endif
   Array<CDT_data> cluster_subdivided(clinfo.tot_cluster());
   for (int c : clinfo.index_range()) {
@@ -2636,19 +2943,11 @@ IMesh trimesh_nary_intersect(const IMesh &tm_in,
   std::cout << "subdivided clusters found, time = "
             << cluster_subdivide_time - subdivided_tris_time << "\n";
 #  endif
-  for (int t : tm_clean->face_index_range()) {
-    int c = clinfo.tri_cluster(t);
-    if (c != NO_INDEX) {
-      BLI_assert(tri_subdivided[t].face_size() == 0);
-      tri_subdivided[t] = extract_subdivided_tri(cluster_subdivided[c], *tm_clean, t, arena);
-    }
-    else if (tri_subdivided[t].face_size() == 0) {
-      tri_subdivided[t] = extract_single_tri(*tm_clean, t);
-    }
-  }
+  calc_cluster_tris(tri_subdivided, *tm_clean, clinfo, cluster_subdivided, arena);
 #  ifdef PERFDEBUG
   double extract_time = PIL_check_seconds_timer();
-  std::cout << "triangles extracted, time = " << extract_time - cluster_subdivide_time << "\n";
+  std::cout << "subdivided cluster tris found, time = " << extract_time - cluster_subdivide_time
+            << "\n";
 #  endif
   IMesh combined = union_tri_subdivides(tri_subdivided);
   if (dbg_level > 1) {