source/blender/blenlib/BLI_offset_indices.hh

@@ -6,13 +6,7 @@
 
 #include <algorithm>
 
-namespace blender {
-namespace index_mask {
-class IndexMask;
-}
-using index_mask::IndexMask;
-}  // namespace blender
-
+#include "BLI_index_mask.hh"
 #include "BLI_index_range.hh"
 #include "BLI_span.hh"
 
@@ -166,6 +160,17 @@ inline OffsetIndices<int> gather_selected_offsets(OffsetIndices<int> src_offsets
 {
   return gather_selected_offsets(src_offsets, selection, 0, dst_offsets);
 }
+OffsetIndices<int> gather_selected_offsets(OffsetIndices<int> src_offsets,
+                                           IndexMaskSegment selection,
+                                           int start_offset,
+                                           MutableSpan<int> dst_offsets);
+inline OffsetIndices<int> gather_selected_offsets(OffsetIndices<int> src_offsets,
+                                                  IndexMaskSegment selection,
+                                                  MutableSpan<int> dst_offsets)
+{
+  return gather_selected_offsets(src_offsets, selection, 0, dst_offsets);
+}
+
 /**
  * Create a map from indexed elements to the source indices, in other words from the larger array
  * to the smaller array.

source/blender/blenlib/BLI_polyfill_2d_beautify.h

@@ -53,6 +53,12 @@ float BLI_polyfill_beautify_quad_rotate_calc_ex(const float v1[2],
 #define BLI_polyfill_beautify_quad_rotate_calc(v1, v2, v3, v4) \
   BLI_polyfill_beautify_quad_rotate_calc_ex(v1, v2, v3, v4, false, NULL)
 
+float BLI_polyfill_edge_calc_rotate_beauty__area(const float v1[3],
+                                                 const float v2[3],
+                                                 const float v3[3],
+                                                 const float v4[3],
+                                                 bool lock_degenerate);
+
 /* avoid realloc's when creating new structures for polyfill ngons */
 #define BLI_POLYFILL_ALLOC_NGON_RESERVE 64
 

source/blender/blenlib/intern/offset_indices.cc

@@ -65,6 +65,20 @@ OffsetIndices<int> gather_selected_offsets(const OffsetIndices<int> src_offsets,
   return OffsetIndices<int>(dst_offsets);
 }
 
+OffsetIndices<int> gather_selected_offsets(const OffsetIndices<int> src_offsets,
+                                           const IndexMaskSegment selection,
+                                           const int start_offset,
+                                           MutableSpan<int> dst_offsets)
+{
+  int offset = start_offset;
+  for (const int64_t i : selection.index_range()) {
+    dst_offsets[i] = offset;
+    offset += src_offsets[selection[i]].size();
+  }
+  dst_offsets.last() = offset;
+  return OffsetIndices<int>(dst_offsets);
+}
+
 void build_reverse_map(OffsetIndices<int> offsets, MutableSpan<int> r_map)
 {
   threading::parallel_for(offsets.index_range(), 1024, [&](const IndexRange range) {

source/blender/blenlib/intern/polyfill_2d_beautify.c

@@ -27,6 +27,7 @@
 
 #include "BLI_heap.h"
 #include "BLI_math_geom.h"
+#include "BLI_math_matrix.h"
 #include "BLI_memarena.h"
 
 #include "BLI_polyfill_2d_beautify.h" /* own include */
@@ -175,6 +176,76 @@ float BLI_polyfill_beautify_quad_rotate_calc_ex(const float v1[2],
   return FLT_MAX;
 }
 
+float BLI_polyfill_edge_calc_rotate_beauty__area(const float v1[3],
+                                                 const float v2[3],
+                                                 const float v3[3],
+                                                 const float v4[3],
+                                                 const bool lock_degenerate)
+{
+  /* not a loop (only to be able to break out) */
+  do {
+    float v1_xy[2], v2_xy[2], v3_xy[2], v4_xy[2];
+
+    /* first get the 2d values */
+    {
+      const float eps = 1e-5f;
+      float no_a[3], no_b[3];
+      float no[3];
+      float axis_mat[3][3];
+      float no_scale;
+      cross_tri_v3(no_a, v2, v3, v4);
+      cross_tri_v3(no_b, v2, v4, v1);
+
+      // printf("%p %p %p %p - %p %p\n", v1, v2, v3, v4, e->l->f, e->l->radial_next->f);
+      BLI_assert((ELEM(v1, v2, v3, v4) == false) && (ELEM(v2, v1, v3, v4) == false) &&
+                 (ELEM(v3, v1, v2, v4) == false) && (ELEM(v4, v1, v2, v3) == false));
+
+      add_v3_v3v3(no, no_a, no_b);
+      if (UNLIKELY((no_scale = normalize_v3(no)) == 0.0f)) {
+        break;
+      }
+
+      axis_dominant_v3_to_m3(axis_mat, no);
+      mul_v2_m3v3(v1_xy, axis_mat, v1);
+      mul_v2_m3v3(v2_xy, axis_mat, v2);
+      mul_v2_m3v3(v3_xy, axis_mat, v3);
+      mul_v2_m3v3(v4_xy, axis_mat, v4);
+
+      /**
+       * Check if input faces are already flipped.
+       * Logic for 'signum_i' addition is:
+       *
+       * Accept:
+       * - (1, 1) or (-1, -1): same side (common case).
+       * - (-1/1, 0): one degenerate, OK since we may rotate into a valid state.
+       *
+       * Ignore:
+       * - (-1, 1): opposite winding, ignore.
+       * - ( 0, 0): both degenerate, ignore.
+       *
+       * \note The cross product is divided by 'no_scale'
+       * so the rotation calculation is scale independent.
+       */
+      if (!(signum_i_ex(cross_tri_v2(v2_xy, v3_xy, v4_xy) / no_scale, eps) +
+            signum_i_ex(cross_tri_v2(v2_xy, v4_xy, v1_xy) / no_scale, eps)))
+      {
+        break;
+      }
+    }
+
+    /**
+     * Important to lock degenerate here,
+     * since the triangle pars will be projected into different 2D spaces.
+     * Allowing to rotate out of a degenerate state can flip the faces
+     * (when performed iteratively).
+     */
+    return BLI_polyfill_beautify_quad_rotate_calc_ex(
+        v1_xy, v2_xy, v3_xy, v4_xy, lock_degenerate, NULL);
+  } while (false);
+
+  return FLT_MAX;
+}
+
 static float polyedge_rotate_beauty_calc(const float (*coords)[2],
                                          const struct HalfEdge *edges,
                                          const struct HalfEdge *e_a,

source/blender/blenloader/intern/versioning_400.cc

@@ -1473,6 +1473,103 @@ static void version_principled_bsdf_coat(bNodeTree *ntree)
       ntree, SH_NODE_BSDF_PRINCIPLED, "Clearcoat Normal", "Coat Normal");
 }
 
+static void remove_triangulate_node_min_size_input(bNodeTree *tree)
+{
+  using namespace blender;
+  Set<bNode *> triangulate_nodes;
+  for (bNode *node : tree->all_nodes()) {
+    if (node->type == GEO_NODE_TRIANGULATE) {
+      triangulate_nodes.add(node);
+    }
+  }
+
+  Map<bNodeSocket *, bNodeLink *> input_links;
+  LISTBASE_FOREACH (bNodeLink *, link, &tree->links) {
+    if (triangulate_nodes.contains(link->tonode)) {
+      input_links.add_new(link->tosock, link);
+    }
+  }
+
+  for (bNode *triangulate : triangulate_nodes) {
+    bNodeSocket *selection = nodeFindSocket(triangulate, SOCK_IN, "Selection");
+    bNodeSocket *min_verts = nodeFindSocket(triangulate, SOCK_IN, "Minimum Vertices");
+    if (!min_verts) {
+      /* Make versioning idempotent. */
+      continue;
+    }
+    const int old_min_verts = static_cast<bNodeSocketValueInt *>(min_verts->default_value)->value;
+    if (!input_links.contains(min_verts) && old_min_verts <= 4) {
+      continue;
+    }
+    bNode *corners_of_face = nodeAddNode(nullptr, tree, "GeometryNodeCornersOfFace");
+    corners_of_face->locx = triangulate->locx - 200;
+    corners_of_face->locy = triangulate->locy - 50;
+    corners_of_face->parent = triangulate->parent;
+    LISTBASE_FOREACH (bNodeSocket *, socket, &corners_of_face->inputs) {
+      socket->flag |= SOCK_HIDDEN;
+    }
+    LISTBASE_FOREACH (bNodeSocket *, socket, &corners_of_face->outputs) {
+      if (!STREQ(socket->identifier, "Total")) {
+        socket->flag |= SOCK_HIDDEN;
+      }
+    }
+
+    bNode *greater_or_equal = nodeAddNode(nullptr, tree, "FunctionNodeCompare");
+    greater_or_equal->locx = triangulate->locx - 100;
+    greater_or_equal->locy = triangulate->locy - 50;
+    greater_or_equal->parent = triangulate->parent;
+    greater_or_equal->flag &= ~NODE_OPTIONS;
+    NodeFunctionCompare *storage = static_cast<NodeFunctionCompare *>(greater_or_equal->storage);
+    storage->operation = NODE_COMPARE_GREATER_EQUAL;
+    storage->data_type = SOCK_INT;
+    nodeAddLink(tree,
+                corners_of_face,
+                nodeFindSocket(corners_of_face, SOCK_OUT, "Total"),
+                greater_or_equal,
+                nodeFindSocket(greater_or_equal, SOCK_IN, "A_INT"));
+    if (bNodeLink **min_verts_link = input_links.lookup_ptr(min_verts)) {
+      (*min_verts_link)->tonode = greater_or_equal;
+      (*min_verts_link)->tosock = nodeFindSocket(greater_or_equal, SOCK_IN, "B_INT");
+    }
+    else {
+      bNodeSocket *new_min_verts = nodeFindSocket(greater_or_equal, SOCK_IN, "B_INT");
+      static_cast<bNodeSocketValueInt *>(new_min_verts->default_value)->value = old_min_verts;
+    }
+
+    if (bNodeLink **selection_link = input_links.lookup_ptr(selection)) {
+      bNode *boolean_and = nodeAddNode(nullptr, tree, "FunctionNodeBooleanMath");
+      boolean_and->locx = triangulate->locx - 75;
+      boolean_and->locy = triangulate->locy - 50;
+      boolean_and->parent = triangulate->parent;
+      boolean_and->flag &= ~NODE_OPTIONS;
+      boolean_and->custom1 = NODE_BOOLEAN_MATH_AND;
+
+      (*selection_link)->tonode = boolean_and;
+      (*selection_link)->tosock = nodeFindSocket(boolean_and, SOCK_IN, "Boolean");
+      nodeAddLink(tree,
+                  greater_or_equal,
+                  nodeFindSocket(greater_or_equal, SOCK_OUT, "Result"),
+                  boolean_and,
+                  nodeFindSocket(boolean_and, SOCK_IN, "Boolean_001"));
+
+      nodeAddLink(tree,
+                  boolean_and,
+                  nodeFindSocket(boolean_and, SOCK_OUT, "Boolean"),
+                  triangulate,
+                  selection);
+    }
+    else {
+      nodeAddLink(tree,
+                  greater_or_equal,
+                  nodeFindSocket(greater_or_equal, SOCK_OUT, "Result"),
+                  triangulate,
+                  selection);
+    }
+
+    /* Make versioning idempotent. */
+    nodeRemoveSocket(tree, triangulate, min_verts);
+  }
+}
 /* Convert specular tint in Principled BSDF. */
 static void version_principled_bsdf_specular_tint(bNodeTree *ntree)
 {
@@ -2575,6 +2672,12 @@ void blo_do_versions_400(FileData *fd, Library * /*lib*/, Main *bmain)
         versioning_nodes_dynamic_sockets_2(*ntree);
       }
     }
+
+    LISTBASE_FOREACH (bNodeTree *, ntree, &bmain->nodetrees) {
+      if (ntree->type == NTREE_GEOMETRY) {
+        remove_triangulate_node_min_size_input(ntree);
+      }
+    }
   }
 
   /**

source/blender/bmesh/tools/bmesh_beautify.cc

@@ -135,76 +135,6 @@ static void erot_state_alternate(const BMEdge *e, EdRotState *e_state)
 /* -------------------------------------------------------------------- */
 /* Calculate the improvement of rotating the edge */
 
-static float bm_edge_calc_rotate_beauty__area(const float v1[3],
-                                              const float v2[3],
-                                              const float v3[3],
-                                              const float v4[3],
-                                              const bool lock_degenerate)
-{
-  /* not a loop (only to be able to break out) */
-  do {
-    float v1_xy[2], v2_xy[2], v3_xy[2], v4_xy[2];
-
-    /* first get the 2d values */
-    {
-      const float eps = 1e-5;
-      float no_a[3], no_b[3];
-      float no[3];
-      float axis_mat[3][3];
-      float no_scale;
-      cross_tri_v3(no_a, v2, v3, v4);
-      cross_tri_v3(no_b, v2, v4, v1);
-
-      // printf("%p %p %p %p - %p %p\n", v1, v2, v3, v4, e->l->f, e->l->radial_next->f);
-      BLI_assert((ELEM(v1, v2, v3, v4) == false) && (ELEM(v2, v1, v3, v4) == false) &&
-                 (ELEM(v3, v1, v2, v4) == false) && (ELEM(v4, v1, v2, v3) == false));
-
-      add_v3_v3v3(no, no_a, no_b);
-      if (UNLIKELY((no_scale = normalize_v3(no)) == 0.0f)) {
-        break;
-      }
-
-      axis_dominant_v3_to_m3(axis_mat, no);
-      mul_v2_m3v3(v1_xy, axis_mat, v1);
-      mul_v2_m3v3(v2_xy, axis_mat, v2);
-      mul_v2_m3v3(v3_xy, axis_mat, v3);
-      mul_v2_m3v3(v4_xy, axis_mat, v4);
-
-      /**
-       * Check if input faces are already flipped.
-       * Logic for 'signum_i' addition is:
-       *
-       * Accept:
-       * - (1, 1) or (-1, -1): same side (common case).
-       * - (-1/1, 0): one degenerate, OK since we may rotate into a valid state.
-       *
-       * Ignore:
-       * - (-1, 1): opposite winding, ignore.
-       * - ( 0, 0): both degenerate, ignore.
-       *
-       * \note The cross product is divided by 'no_scale'
-       * so the rotation calculation is scale independent.
-       */
-      if (!(signum_i_ex(cross_tri_v2(v2_xy, v3_xy, v4_xy) / no_scale, eps) +
-            signum_i_ex(cross_tri_v2(v2_xy, v4_xy, v1_xy) / no_scale, eps)))
-      {
-        break;
-      }
-    }
-
-    /**
-     * Important to lock degenerate here,
-     * since the triangle pars will be projected into different 2D spaces.
-     * Allowing to rotate out of a degenerate state can flip the faces
-     * (when performed iteratively).
-     */
-    return BLI_polyfill_beautify_quad_rotate_calc_ex(
-        v1_xy, v2_xy, v3_xy, v4_xy, lock_degenerate, nullptr);
-  } while (false);
-
-  return FLT_MAX;
-}
-
 static float bm_edge_calc_rotate_beauty__angle(const float v1[3],
                                                const float v2[3],
                                                const float v3[3],
@@ -256,7 +186,7 @@ float BM_verts_calc_rotate_beauty(const BMVert *v1,
 
     switch (method) {
       case 0:
-        return bm_edge_calc_rotate_beauty__area(
+        return BLI_polyfill_edge_calc_rotate_beauty__area(
             v1->co, v2->co, v3->co, v4->co, flag & EDGE_RESTRICT_DEGENERATE);
       default:
         return bm_edge_calc_rotate_beauty__angle(v1->co, v2->co, v3->co, v4->co);

source/blender/geometry/CMakeLists.txt

@@ -30,6 +30,7 @@ set(SRC
   intern/mesh_split_edges.cc
   intern/mesh_to_curve_convert.cc
   intern/mesh_to_volume.cc
+  intern/mesh_triangulate.cc
   intern/point_merge_by_distance.cc
   intern/points_to_volume.cc
   intern/randomize.cc
@@ -57,6 +58,7 @@ set(SRC
   GEO_mesh_split_edges.hh
   GEO_mesh_to_curve.hh
   GEO_mesh_to_volume.hh
+  GEO_mesh_triangulate.hh
   GEO_point_merge_by_distance.hh
   GEO_points_to_volume.hh
   GEO_randomize.hh
@@ -74,6 +76,7 @@ set(LIB
   bf_blenkernel
   PRIVATE bf::blenlib
   PRIVATE bf::dna
+  PRIVATE bf::intern::atomic
   PRIVATE bf::intern::guardedalloc
 )
 

source/blender/geometry/GEO_mesh_triangulate.hh

@@ -0,0 +1,49 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include <optional>
+
+#include "BLI_index_mask.hh"
+
+#include "BKE_attribute.hh"
+
+struct Mesh;
+
+namespace blender::geometry {
+
+/** \warning Values are saved in files. */
+enum class TriangulateNGonMode {
+  /** Add a "beauty" pass on top of the standard ear-clipping algorithm. */
+  Beauty = 0,
+  EarClip = 1,
+};
+
+/** \warning Values are saved in files. */
+enum class TriangulateQuadMode {
+  /** Complex method to determine the best looking edge. */
+  Beauty = 0,
+  /** Create a new edge from the first corner to the last. */
+  Fixed = 1,
+  /** Create a new edge from the second corner to the third. */
+  Alternate = 2,
+  /** Create a new edge along the shortest diagonal. */
+  ShortEdge = 3,
+  /** Create a new edge along the longest diagonal. */
+  LongEdge = 4,
+};
+
+/**
+ * \return #std::nullopt if the mesh is not changed (when every selected face is already a
+ * triangle).
+ */
+std::optional<Mesh *> mesh_triangulate(
+    const Mesh &src_mesh,
+    const IndexMask &selection,
+    TriangulateNGonMode ngon_mode,
+    TriangulateQuadMode quad_mode,
+    const bke::AnonymousAttributePropagationInfo &propagation_info);
+
+}  // namespace blender::geometry

source/blender/geometry/intern/mesh_copy_selection.cc

@@ -383,6 +383,8 @@ std::optional<Mesh *> mesh_copy_selection_keep_verts(
       [&]() {
         bke::copy_attributes(
             src_attributes, bke::AttrDomain::Point, propagation_info, {}, dst_attributes);
+        CustomData_merge(
+            &src_mesh.vert_data, &dst_mesh->vert_data, CD_MASK_ORIGINDEX, src_mesh.verts_num);
         bke::gather_attributes(src_attributes,
                                bke::AttrDomain::Edge,
                                propagation_info,
@@ -464,8 +466,12 @@ std::optional<Mesh *> mesh_copy_selection_keep_edges(
 
   bke::copy_attributes(
       src_attributes, bke::AttrDomain::Point, propagation_info, {}, dst_attributes);
+  CustomData_merge(
+      &src_mesh.vert_data, &dst_mesh->vert_data, CD_MASK_ORIGINDEX, src_mesh.verts_num);
   bke::copy_attributes(
       src_attributes, bke::AttrDomain::Edge, propagation_info, {}, dst_attributes);
+  CustomData_merge(
+      &src_mesh.edge_data, &dst_mesh->edge_data, CD_MASK_ORIGINDEX, src_mesh.edges_num);
   bke::gather_attributes(
       src_attributes, bke::AttrDomain::Face, propagation_info, {}, face_mask, dst_attributes);
   bke::gather_attributes_group_to_group(src_attributes,

source/blender/geometry/intern/mesh_triangulate.cc

@@ -0,0 +1,890 @@
+/* SPDX-FileCopyrightText: 2023 Blender Authors
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "atomic_ops.h"
+
+#include "BLI_array_utils.hh"
+#include "BLI_enumerable_thread_specific.hh"
+#include "BLI_index_mask.hh"
+#include "BLI_math_geom.h"
+#include "BLI_math_matrix.h"
+#include "BLI_ordered_edge.hh"
+#include "BLI_polyfill_2d.h"
+#include "BLI_polyfill_2d_beautify.h"
+#include "BLI_vector_set.hh"
+
+#include "BLI_heap.h"
+#include "BLI_memarena.h"
+
+#include "BKE_attribute.hh"
+#include "BKE_attribute_math.hh"
+#include "BKE_customdata.hh"
+#include "BKE_mesh.hh"
+#include "BKE_mesh_mapping.hh"
+
+#include "GEO_mesh_triangulate.hh"
+
+namespace blender::geometry {
+
+static void gather(const Span<int> src, const Span<int16_t> indices, MutableSpan<int> dst)
+{
+  for (const int i : indices.index_range()) {
+    dst[i] = src[indices[i]];
+  }
+}
+
+static Span<int> gather_or_reference(const Span<int> src,
+                                     const Span<int16_t> indices,
+                                     Vector<int> &dst)
+{
+  if (unique_sorted_indices::non_empty_is_range(indices)) {
+    return src.slice(indices[0], indices.size());
+  }
+  dst.reinitialize(indices.size());
+  gather(src, indices, dst);
+  return dst.as_span();
+}
+
+static Span<int> gather_or_reference(const Span<int> src,
+                                     const IndexMaskSegment mask,
+                                     Vector<int> &dst)
+{
+  return gather_or_reference(src.drop_front(mask.offset()), mask.base_span(), dst);
+}
+
+/**
+ * If a significant number of Ngons are selected (> 25% of the faces), then use the
+ * face normals cache, in case the cache is persistent (or already calculated).
+ */
+static Span<float3> face_normals_if_worthwhile(const Mesh &src_mesh, const int selection_size)
+{
+  if (src_mesh.runtime->face_normals_cache.is_cached()) {
+    return src_mesh.face_normals();
+  }
+  if (selection_size > src_mesh.faces_num / 4) {
+    return src_mesh.face_normals();
+  }
+  return {};
+}
+
+static void copy_loose_vert_hint(const Mesh &src, Mesh &dst)
+{
+  const auto &src_cache = src.runtime->loose_verts_cache;
+  if (src_cache.is_cached() && src_cache.data().count == 0) {
+    dst.tag_loose_verts_none();
+  }
+}
+
+static void copy_loose_edge_hint(const Mesh &src, Mesh &dst)
+{
+  const auto &src_cache = src.runtime->loose_edges_cache;
+  if (src_cache.is_cached() && src_cache.data().count == 0) {
+    dst.tag_loose_edges_none();
+  }
+}
+
+static Mesh *create_mesh_no_attributes(const Mesh &params_mesh,
+                                       const int verts_num,
+                                       const int edges_num,
+                                       const int faces_num,
+                                       const int corners_num)
+{
+  Mesh *mesh = BKE_mesh_new_nomain(0, 0, faces_num, 0);
+  BKE_mesh_copy_parameters_for_eval(mesh, &params_mesh);
+  CustomData_free_layer_named(&mesh->vert_data, "position", 0);
+  CustomData_free_layer_named(&mesh->edge_data, ".edge_verts", 0);
+  CustomData_free_layer_named(&mesh->corner_data, ".corner_vert", 0);
+  CustomData_free_layer_named(&mesh->corner_data, ".corner_edge", 0);
+  mesh->verts_num = verts_num;
+  mesh->edges_num = edges_num;
+  mesh->corners_num = corners_num;
+  return mesh;
+}
+
+static OffsetIndices<int> calc_faces(const OffsetIndices<int> src_faces,
+                                     const IndexMask &unselected,
+                                     MutableSpan<int> offsets)
+{
+  MutableSpan<int> new_tri_offsets = offsets.drop_back(unselected.size());
+  offset_indices::fill_constant_group_size(3, new_tri_offsets.first(), new_tri_offsets);
+  offset_indices::gather_selected_offsets(
+      src_faces, unselected, new_tri_offsets.last(), offsets.take_back(unselected.size() + 1));
+  return OffsetIndices<int>(offsets);
+}
+
+namespace quad {
+
+/**
+ *  #Edge_0_2       #Edge_1_3
+ * 3 ------- 2     3 ------- 2
+ * | 1     / |     | \     1 |
+ * |     /   |     |   \     |
+ * |   /     |     |     \   |
+ * | /     0 |     | 0     \ |
+ * 0 ------- 1     0 ------- 1
+ */
+enum class QuadDirection : int8_t {
+  Edge_0_2 = 0,
+  Edge_1_3 = 1,
+};
+
+/**
+ * \note This behavior is meant to be the same as #BM_verts_calc_rotate_beauty.
+ * The order of vertices requires special attention.
+ */
+static QuadDirection calc_quad_direction_beauty(const float3 &v0,
+                                                const float3 &v1,
+                                                const float3 &v2,
+                                                const float3 &v3)
+{
+  const int flip_flag = is_quad_flip_v3(v1, v2, v3, v0);
+  if (UNLIKELY(flip_flag & (1 << 0))) {
+    return QuadDirection::Edge_0_2;
+  }
+  if (UNLIKELY(flip_flag & (1 << 1))) {
+    return QuadDirection::Edge_1_3;
+  }
+  return BLI_polyfill_edge_calc_rotate_beauty__area(v1, v2, v3, v0, false) > 0.0f ?
+             QuadDirection::Edge_0_2 :
+             QuadDirection::Edge_1_3;
+}
+
+static void calc_quad_directions(const Span<float3> positions,
+                                 const Span<int> face_offsets,
+                                 const Span<int> corner_verts,
+                                 const TriangulateQuadMode quad_mode,
+                                 MutableSpan<QuadDirection> directions)
+{
+  if (quad_mode == TriangulateQuadMode::Fixed) {
+    directions.fill(QuadDirection::Edge_0_2);
+  }
+  else if (quad_mode == TriangulateQuadMode::Alternate) {
+    directions.fill(QuadDirection::Edge_1_3);
+  }
+  else if (quad_mode == TriangulateQuadMode::Beauty) {
+    for (const int i : face_offsets.index_range()) {
+      const Span<int> verts = corner_verts.slice(face_offsets[i], 4);
+      directions[i] = calc_quad_direction_beauty(
+          positions[verts[0]], positions[verts[1]], positions[verts[2]], positions[verts[3]]);
+    }
+  }
+  else {
+    const QuadDirection long_dir = quad_mode == TriangulateQuadMode::ShortEdge ?
+                                       QuadDirection::Edge_0_2 :
+                                       QuadDirection::Edge_1_3;
+    const QuadDirection short_dir = quad_mode == TriangulateQuadMode::ShortEdge ?
+                                        QuadDirection::Edge_1_3 :
+                                        QuadDirection::Edge_0_2;
+    for (const int i : face_offsets.index_range()) {
+      const Span<int> verts = corner_verts.slice(face_offsets[i], 4);
+      const float dist_0_2 = math::distance_squared(positions[verts[0]], positions[verts[2]]);
+      const float dist_1_3 = math::distance_squared(positions[verts[1]], positions[verts[3]]);
+      directions[i] = dist_0_2 < dist_1_3 ? long_dir : short_dir;
+    }
+  }
+}
+
+static void calc_corner_maps(const Span<int> face_offsets,
+                             const Span<QuadDirection> directions,
+                             MutableSpan<int> corner_map)
+{
+  for (const int i : face_offsets.index_range()) {
+    MutableSpan<int> quad_map = corner_map.slice(6 * i, 6);
+    /* These corner orders give new edges based on the first vertex of each triangle. */
+    switch (directions[i]) {
+      case QuadDirection::Edge_0_2:
+        quad_map.copy_from({2, 0, 1, 0, 2, 3});
+        break;
+      case QuadDirection::Edge_1_3:
+        quad_map.copy_from({1, 3, 0, 3, 1, 2});
+        break;
+    }
+    const int face_start = face_offsets[i];
+    for (int &i : quad_map) {
+      i += face_start;
+    }
+  }
+}
+
+static void calc_corner_maps(const Span<float3> positions,
+                             const OffsetIndices<int> src_faces,
+                             const Span<int> src_corner_verts,
+                             const IndexMask &quads,
+                             const TriangulateQuadMode quad_mode,
+                             MutableSpan<int> corner_map)
+{
+  struct TLS {
+    Vector<int> offsets;
+    Vector<QuadDirection> directions;
+  };
+  threading::EnumerableThreadSpecific<TLS> tls;
+
+  quads.foreach_segment(GrainSize(1024), [&](const IndexMaskSegment quads, const int64_t pos) {
+    TLS &data = tls.local();
+    data.directions.reinitialize(quads.size());
+
+    /* Find the offsets of each face in the local selection. We can gather them together even if
+     * they aren't contiguous because we only need to know the start of each face; the size is
+     * just 4. */
+    const Span<int> offsets = gather_or_reference(src_faces.data(), quads, data.offsets);
+    calc_quad_directions(positions, offsets, src_corner_verts, quad_mode, data.directions);
+    const IndexRange corners(pos * 6, offsets.size() * 6);
+    quad::calc_corner_maps(offsets, data.directions, corner_map.slice(corners));
+  });
+}
+
+/**
+ * Each triangulated quad creates one additional edge in the result mesh, between the two
+ * triangles. The corner_verts are just the corners of the quads, and the edges are just the new
+ * edges for these quads.
+ */
+static void calc_edges(const Span<int> quad_corner_verts, MutableSpan<int2> new_quad_edges)
+{
+  const int quads_num = quad_corner_verts.size() / 6;
+  for (const int i : IndexRange(quads_num)) {
+    const Span<int> verts = quad_corner_verts.slice(6 * i, 6);
+    /* Use the first vertex of each triangle. */
+    new_quad_edges[i] = int2(verts[0], verts[1]);
+  }
+}
+
+static void calc_quad_corner_edges(const Span<int> src_corner_edges,
+                                   const Span<int> corner_map,
+                                   const int edges_start,
+                                   MutableSpan<int> corner_edges)
+{
+  /* Each triangle starts at the new edge and winds in the same order as corner vertices
+   * described by the corner map. */
+  for (const int tri : IndexRange(corner_map.size() / 3)) {
+    const int tri_start = 3 * tri;
+    corner_edges[tri_start] = edges_start + tri / 2;
+    corner_edges[tri_start + 1] = src_corner_edges[corner_map[tri_start + 1]];
+    corner_edges[tri_start + 2] = src_corner_edges[corner_map[tri_start + 2]];
+  }
+}
+
+static void calc_edges(const Span<int> src_corner_edges,
+                       const Span<int> corner_map,
+                       const Span<int> corner_verts,
+                       const int edges_start,
+                       MutableSpan<int2> edges,
+                       MutableSpan<int> quad_corner_edges)
+{
+  const int quads_num = corner_map.size() / 6;
+  threading::parallel_for(IndexRange(quads_num), 1024, [&](const IndexRange quads) {
+    const IndexRange corners(quads.start() * 6, quads.size() * 6);
+    calc_edges(corner_verts.slice(corners), edges.slice(quads));
+    calc_quad_corner_edges(src_corner_edges,
+                           corner_map.slice(corners),
+                           edges_start + quads.start(),
+                           quad_corner_edges.slice(corners));
+  });
+}
+
+template<typename T>
+static void copy_face_data_to_tris(const Span<T> src, const IndexMask &quads, MutableSpan<T> dst)
+{
+  quads.foreach_index_optimized<int>([&](const int src_i, const int dst_i) {
+    dst[2 * dst_i + 0] = src[src_i];
+    dst[2 * dst_i + 1] = src[src_i];
+  });
+}
+
+static void copy_face_data_to_tris(const GSpan src, const IndexMask &quads, GMutableSpan dst)
+{
+  bke::attribute_math::convert_to_static_type(src.type(), [&](auto dummy) {
+    using T = decltype(dummy);
+    copy_face_data_to_tris(src.typed<T>(), quads, dst.typed<T>());
+  });
+}
+
+}  // namespace quad
+
+namespace ngon {
+
+static OffsetIndices<int> calc_tris_by_ngon(const OffsetIndices<int> src_faces,
+                                            const IndexMask &ngons,
+                                            MutableSpan<int> face_offset_data)
+{
+  ngons.foreach_index(GrainSize(2048), [&](const int face, const int mask) {
+    face_offset_data[mask] = bke::mesh::face_triangles_num(src_faces[face].size());
+  });
+  return offset_indices::accumulate_counts_to_offsets(face_offset_data);
+}
+
+static OffsetIndices<int> calc_edges_by_ngon(const OffsetIndices<int> src_faces,
+                                             const IndexMask &selection,
+                                             MutableSpan<int> edge_offset_data)
+{
+  selection.foreach_index(GrainSize(2048), [&](const int face, const int mask) {
+    /* The number of new inner edges for each face is the number of corners - 3. */
+    edge_offset_data[mask] = src_faces[face].size() - 3;
+  });
+  return offset_indices::accumulate_counts_to_offsets(edge_offset_data);
+}
+
+static void calc_corner_maps(const Span<float3> positions,
+                             const OffsetIndices<int> src_faces,
+                             const Span<int> src_corner_verts,
+                             const Span<float3> face_normals,
+                             const IndexMask &ngons,
+                             const OffsetIndices<int> tris_by_ngon,
+                             const TriangulateNGonMode ngon_mode,
+                             MutableSpan<int> corner_map)
+{
+  struct TLS {
+    Vector<float3x3> projections;
+    Array<int> offset_data;
+    Vector<float2> projected_positions;
+
+    /* Only used for the "Beauty" method. */
+    MemArena *arena = nullptr;
+    Heap *heap = nullptr;
+
+    ~TLS()
+    {
+      if (arena) {
+        BLI_memarena_free(arena);
+      }
+      if (heap) {
+        BLI_heap_free(heap, nullptr);
+      }
+    }
+  };
+  threading::EnumerableThreadSpecific<TLS> tls;
+
+  ngons.foreach_segment(GrainSize(128), [&](const IndexMaskSegment ngons, const int pos) {
+    TLS &data = tls.local();
+
+    /* In order to simplify and "parallelize" the next loops, gather offsets used to group an array
+     * large enough for all the local face corners. */
+    data.offset_data.reinitialize(ngons.size() + 1);
+    const OffsetIndices local_corner_offsets = offset_indices::gather_selected_offsets(
+        src_faces, ngons, data.offset_data);
+
+    /* Use face normals to build projection matrices to make the face positions 2D. */
+    data.projections.reinitialize(ngons.size());
+    MutableSpan<float3x3> projections = data.projections;
+    if (face_normals.is_empty()) {
+      for (const int i : ngons.index_range()) {
+        const IndexRange src_face = src_faces[ngons[i]];
+        const Span<int> face_verts = src_corner_verts.slice(src_face);
+        const float3 normal = bke::mesh::face_normal_calc(positions, face_verts);
+        axis_dominant_v3_to_m3_negate(projections[i].ptr(), normal);
+      }
+    }
+    else {
+      for (const int i : ngons.index_range()) {
+        axis_dominant_v3_to_m3_negate(projections[i].ptr(), face_normals[ngons[i]]);
+      }
+    }
+
+    /* Project the face positions into 2D using the matrices calculated above. */
+    data.projected_positions.reinitialize(local_corner_offsets.total_size());
+    MutableSpan<float2> projected_positions = data.projected_positions;
+    for (const int i : ngons.index_range()) {
+      const IndexRange src_face = src_faces[ngons[i]];
+      const Span<int> face_verts = src_corner_verts.slice(src_face);
+      const float3x3 &matrix = projections[i];
+
+      MutableSpan<float2> positions_2d = projected_positions.slice(local_corner_offsets[i]);
+      for (const int i : face_verts.index_range()) {
+        mul_v2_m3v3(positions_2d[i], matrix.ptr(), positions[face_verts[i]]);
+      }
+    }
+
+    if (ngon_mode == TriangulateNGonMode::Beauty) {
+      if (!data.arena) {
+        data.arena = BLI_memarena_new(BLI_POLYFILL_ARENA_SIZE, __func__);
+      }
+      if (!data.heap) {
+        data.heap = BLI_heap_new_ex(BLI_POLYFILL_ALLOC_NGON_RESERVE);
+      }
+    }
+
+    /* Calculate the triangulation of corners indices local to each face. */
+    for (const int i : ngons.index_range()) {
+      const Span<float2> positions_2d = projected_positions.slice(local_corner_offsets[i]);
+      const IndexRange tris_range = tris_by_ngon[pos + i];
+      MutableSpan<int> map = corner_map.slice(tris_range.start() * 3, tris_range.size() * 3);
+      BLI_polyfill_calc(reinterpret_cast<const float(*)[2]>(positions_2d.data()),
+                        positions_2d.size(),
+                        1,
+                        reinterpret_cast<uint(*)[3]>(map.data()));
+      if (ngon_mode == TriangulateNGonMode::Beauty) {
+        BLI_polyfill_beautify(reinterpret_cast<const float(*)[2]>(positions_2d.data()),
+                              positions_2d.size(),
+                              reinterpret_cast<uint(*)[3]>(map.data()),
+                              data.arena,
+                              data.heap);
+        BLI_memarena_clear(data.arena);
+      }
+    }
+
+    /* "Globalize" the triangulation created above so the map source indices reference _all_ of the
+     * source vertices, not just within the source face. */
+    for (const int i : ngons.index_range()) {
+      const IndexRange src_face = src_faces[ngons[i]];
+      const IndexRange tris_range = tris_by_ngon[pos + i];
+      MutableSpan<int> map = corner_map.slice(tris_range.start() * 3, tris_range.size() * 3);
+      for (int &vert : map) {
+        vert += src_face.start();
+      }
+    }
+  });
+}
+
+static void calc_inner_tri_edges(const IndexRange src_face,
+                                 const Span<int> src_corner_verts,
+                                 const Span<int> src_corner_edges,
+                                 const Span<int3> corner_tris,
+                                 const int edges_start,
+                                 MutableSpan<int> corner_edges,
+                                 VectorSet<OrderedEdge> &deduplication)
+{
+  const OrderedEdge last_edge(int(src_face.first()), int(src_face.last()));
+  auto add_edge = [&](const OrderedEdge corner_edge) -> int {
+    if (corner_edge == last_edge) {
+      return src_corner_edges[src_face.last()];
+    }
+    if (corner_edge.v_high == corner_edge.v_low + 1) {
+      return src_corner_edges[corner_edge.v_low];
+    }
+    const OrderedEdge vert_edge(src_corner_verts[corner_edge.v_low],
+                                src_corner_verts[corner_edge.v_high]);
+    return edges_start + deduplication.index_of_or_add(vert_edge);
+  };
+
+  for (const int i : corner_tris.index_range()) {
+    const int3 tri = corner_tris[i];
+    corner_edges[3 * i + 0] = add_edge({tri[0], tri[1]});
+    corner_edges[3 * i + 1] = add_edge({tri[1], tri[2]});
+    corner_edges[3 * i + 2] = add_edge({tri[2], tri[0]});
+  }
+}
+
+static void calc_edges(const OffsetIndices<int> src_faces,
+                       const Span<int> src_corner_verts,
+                       const Span<int> src_corner_edges,
+                       const IndexMask &ngons,
+                       const OffsetIndices<int> tris_by_ngon,
+                       const OffsetIndices<int> edges_by_ngon,
+                       const OffsetIndices<int> faces,
+                       const int edges_start,
+                       const Span<int> corner_verts,
+                       MutableSpan<int2> inner_edges,
+                       MutableSpan<int> corner_edges)
+{
+  threading::EnumerableThreadSpecific<VectorSet<OrderedEdge>> tls;
+  ngons.foreach_segment(GrainSize(128), [&](const IndexMaskSegment ngons, const int pos) {
+    VectorSet<OrderedEdge> &deduplication = tls.local();
+    for (const int i : ngons.index_range()) {
+      const IndexRange edges = edges_by_ngon[pos + i];
+      const IndexRange corners = faces[tris_by_ngon[pos + i]];
+      deduplication.clear();
+      calc_inner_tri_edges(src_faces[ngons[i]],
+                           src_corner_verts,
+                           src_corner_edges,
+                           corner_verts.slice(corners).cast<int3>(),
+                           edges_start + edges.start(),
+                           corner_edges.slice(corners),
+                           deduplication);
+      inner_edges.slice(edges).copy_from(deduplication.as_span().cast<int2>());
+    }
+  });
+}
+
+}  // namespace ngon
+
+namespace deduplication {
+
+static GroupedSpan<int> build_vert_to_tri_map(const int verts_num,
+                                              const Span<int3> vert_tris,
+                                              Array<int> &r_offsets,
+                                              Array<int> &r_indices)
+{
+  r_offsets = Array<int>(verts_num + 1, 0);
+  offset_indices::build_reverse_offsets(vert_tris.cast<int>(), r_offsets);
+  const OffsetIndices offsets(r_offsets.as_span());
+
+  r_indices.reinitialize(offsets.total_size());
+  int *counts = MEM_cnew_array<int>(size_t(offsets.size()), __func__);
+  BLI_SCOPED_DEFER([&]() { MEM_freeN(counts); })
+  threading::parallel_for(vert_tris.index_range(), 1024, [&](const IndexRange range) {
+    for (const int tri : range) {
+      for (const int vert : {vert_tris[tri][0], vert_tris[tri][1], vert_tris[tri][2]}) {
+        const int index_in_group = atomic_fetch_and_add_int32(&counts[vert], 1);
+        r_indices[offsets[vert][index_in_group]] = tri;
+      }
+    }
+  });
+
+  return {r_offsets.as_span(), r_indices.as_span()};
+}
+
+/**
+ * To avoid adding duplicate faces to the mesh without complicating the triangulation code to
+ * support that unlikely case, check if triangles (which are all unselected) have an equivalent
+ * newly created triangle, and don't copy them to the result mesh if so.
+ */
+static IndexMask calc_unselected_faces(const Mesh &mesh,
+                                       const OffsetIndices<int> src_faces,
+                                       const Span<int> src_corner_verts,
+                                       const IndexMask &selection,
+                                       const Span<int> corner_map,
+                                       IndexMaskMemory &memory)
+{
+  const IndexMask unselected = selection.complement(src_faces.index_range(), memory);
+  if (mesh.no_overlapping_topology()) {
+    return unselected;
+  }
+  const bool has_tris = threading::parallel_reduce(
+      unselected.index_range(),
+      4096,
+      false,
+      [&](const IndexRange range, const bool init) {
+        if (init) {
+          return init;
+        }
+        const IndexMask sliced_mask = unselected.slice(range);
+        for (const int64_t segment_i : IndexRange(sliced_mask.segments_num())) {
+          const IndexMaskSegment segment = sliced_mask.segment(segment_i);
+          for (const int i : segment) {
+            if (src_faces[i].size() == 3) {
+              return true;
+            }
+          }
+        }
+        return false;
+      },
+      std::logical_or());
+  if (!has_tris) {
+    return unselected;
+  }
+  Array<int3> vert_tris(corner_map.size() / 3);
+  bke::attribute_math::gather(
+      src_corner_verts, corner_map, vert_tris.as_mutable_span().cast<int>());
+
+  Array<int> vert_to_tri_offsets;
+  Array<int> vert_to_tri_indices;
+  const GroupedSpan<int> vert_to_tri = build_vert_to_tri_map(
+      mesh.verts_num, corner_map.cast<int3>(), vert_to_tri_offsets, vert_to_tri_indices);
+  return IndexMask::from_predicate(unselected, GrainSize(1024), memory, [&](const int i) {
+    const Span<int> face_verts = src_corner_verts.slice(src_faces[i]);
+    if (face_verts.size() != 3) {
+      return true;
+    }
+    /* TODO: Sorting the three values with a few comparisons would be faster than a #Set. */
+    const Set<int, 3> vert_set(face_verts);
+    return std::none_of(face_verts.begin(), face_verts.end(), [&](const int vert) {
+      return std::none_of(vert_to_tri[vert].begin(), vert_to_tri[vert].end(), [&](const int tri) {
+        return Set<int, 3>(Span(&vert_tris[tri].x, 3)) == vert_set;
+      });
+    });
+  });
+}
+
+static std::optional<int> find_edge_duplicate(const GroupedSpan<int> vert_to_edge_map,
+                                              const Span<int2> edges,
+                                              const OrderedEdge edge)
+{
+  for (const int vert : {edge.v_low, edge.v_high}) {
+    for (const int src_edge : vert_to_edge_map[vert]) {
+      if (OrderedEdge(edges[src_edge]) == edge) {
+        return src_edge;
+      }
+    }
+  }
+  return std::nullopt;
+}
+
+/**
+ * Given all the edges on the new mesh, find new edges that are duplicates of existing edges.
+ * If there are any, remove them and references to them in the corner edge array.
+ *
+ * \return The final number of edges in the mesh.
+ */
+static int calc_new_edges(const Mesh &src_mesh,
+                          const Span<int2> src_edges,
+                          const GroupedSpan<int> vert_to_edge,
+                          const IndexRange new_edges_range,
+                          MutableSpan<int2> edges,
+                          MutableSpan<int> corner_edges)
+{
+  if (src_mesh.no_overlapping_topology()) {
+    return edges.size();
+  }
+  const Span<int2> new_edges = edges.slice(new_edges_range);
+  Array<int> duplicate_remap(new_edges.size());
+  threading::parallel_for(new_edges.index_range(), 1024, [&](const IndexRange range) {
+    for (const int i : range) {
+      duplicate_remap[i] = find_edge_duplicate(vert_to_edge, src_edges, new_edges[i]).value_or(-1);
+    }
+  });
+  IndexMaskMemory memory;
+  const IndexMask non_duplicate_new_edges = IndexMask::from_predicate(
+      new_edges.index_range(), GrainSize(4096), memory, [&](const int i) {
+        return duplicate_remap[i] == -1;
+      });
+  if (non_duplicate_new_edges.size() == new_edges.size()) {
+    return edges.size();
+  }
+
+  non_duplicate_new_edges.foreach_index_optimized<int>(
+      GrainSize(4096), [&](const int index, const int pos) {
+        duplicate_remap[index] = pos + new_edges_range.start();
+      });
+  threading::parallel_for(corner_edges.index_range(), 4096, [&](const IndexRange range) {
+    for (const int corner : range) {
+      const int edge = corner_edges[corner];
+      if (edge < new_edges_range.start()) {
+        continue;
+      }
+      const int remap_index = edge - new_edges_range.start();
+      corner_edges[corner] = duplicate_remap[remap_index];
+    }
+  });
+
+  Array<int2> edges_with_duplicates(new_edges);
+  array_utils::gather(edges_with_duplicates.as_span(),
+                      non_duplicate_new_edges,
+                      edges.slice(new_edges_range.start(), non_duplicate_new_edges.size()));
+  return src_edges.size() + non_duplicate_new_edges.size();
+}
+
+}  // namespace deduplication
+
+std::optional<Mesh *> mesh_triangulate(
+    const Mesh &src_mesh,
+    const IndexMask &selection_with_tris,
+    const TriangulateNGonMode ngon_mode,
+    const TriangulateQuadMode quad_mode,
+    const bke::AnonymousAttributePropagationInfo &propagation_info)
+{
+  const Span<float3> positions = src_mesh.vert_positions();
+  const Span<int2> src_edges = src_mesh.edges();
+  const OffsetIndices src_faces = src_mesh.faces();
+  const Span<int> src_corner_verts = src_mesh.corner_verts();
+  const Span<int> src_corner_edges = src_mesh.corner_edges();
+  const bke::AttributeAccessor src_attributes = src_mesh.attributes();
+
+  /* Divide the input selection into separate selections for each face type. This isn't necessary
+   * for correctness, but considering groups of each face type separately simplifies optimizing
+   * for each type. For example, quad triangulation is much simpler than Ngon triangulation. */
+  IndexMaskMemory memory;
+  const IndexMask selection = IndexMask::from_predicate(
+      selection_with_tris, GrainSize(4096), memory, [&](const int i) {
+        return src_faces[i].size() > 3;
+      });
+  const IndexMask quads = IndexMask::from_predicate(
+      selection, GrainSize(4096), memory, [&](const int i) { return src_faces[i].size() == 4; });
+  const IndexMask ngons = IndexMask::from_predicate(
+      selection, GrainSize(4096), memory, [&](const int i) { return src_faces[i].size() > 4; });
+  if (quads.is_empty() && ngons.is_empty()) {
+    /* All selected faces are already triangles. */
+    return std::nullopt;
+  }
+
+  /* Calculate group of triangle indices for each selected Ngon to facilitate calculating them in
+   * parallel later. */
+  Array<int> tri_offsets(ngons.size() + 1);
+  const OffsetIndices tris_by_ngon = ngon::calc_tris_by_ngon(src_faces, ngons, tri_offsets);
+  const int ngon_tris_num = tris_by_ngon.total_size();
+  const int quad_tris_num = quads.size() * 2;
+  const IndexRange tris_range(ngon_tris_num + quad_tris_num);
+  const IndexRange ngon_tris_range = tris_range.take_front(ngon_tris_num);
+  const IndexRange quad_tris_range = tris_range.take_front(quad_tris_num);
+
+  const int ngon_corners_num = tris_by_ngon.total_size() * 3;
+  const int quad_corners_num = quads.size() * 6;
+  const IndexRange tri_corners_range(quad_corners_num + ngon_corners_num);
+  const IndexRange ngon_corners_range = tri_corners_range.take_front(ngon_corners_num);
+  const IndexRange quad_corners_range = tri_corners_range.take_back(quad_corners_num);
+
+  /* Calculate groups of new inner edges for each selected Ngon so they can be filled in parallel
+   * later. */
+  Array<int> edge_offset_data(ngons.size() + 1);
+  const OffsetIndices edges_by_ngon = ngon::calc_edges_by_ngon(src_faces, ngons, edge_offset_data);
+  const int ngon_edges_num = edges_by_ngon.total_size();
+  const int quad_edges_num = quads.size();
+  const IndexRange src_edges_range(0, src_edges.size());
+  const IndexRange tri_edges_range(src_edges_range.one_after_last(),
+                                   ngon_edges_num + quad_edges_num);
+  const IndexRange ngon_edges_range = tri_edges_range.take_front(ngon_edges_num);
+  const IndexRange quad_edges_range = tri_edges_range.take_front(quad_edges_num);
+
+  /* An index map that maps from newly created corners in `tri_corners_range` to original corner
+   * indices. This is used to interpolate `corner_vert` indices and face corner attributes. If
+   * there are no face corner attributes, theoretically the map could be skipped and corner
+   * vertex indices could be interpolated immediately, but that isn't done for simplicity. */
+  Array<int> corner_map(tri_corners_range.size());
+
+  if (!ngons.is_empty()) {
+    ngon::calc_corner_maps(positions,
+                           src_faces,
+                           src_corner_verts,
+                           face_normals_if_worthwhile(src_mesh, ngons.size()),
+                           ngons,
+                           tris_by_ngon,
+                           ngon_mode,
+                           corner_map.as_mutable_span().slice(ngon_corners_range));
+  }
+  if (!quads.is_empty()) {
+    quad::calc_corner_maps(positions,
+                           src_faces,
+                           src_corner_verts,
+                           quads,
+                           quad_mode,
+                           corner_map.as_mutable_span().slice(quad_corners_range));
+  }
+
+  const IndexMask unselected = deduplication::calc_unselected_faces(
+      src_mesh, src_faces, src_corner_verts, selection, corner_map, memory);
+  const IndexRange unselected_range(tris_range.one_after_last(), unselected.size());
+
+  /* Create a mesh with no face corners.
+   *  - We haven't yet counted the number of corners from unselected faces. Creating the final face
+   *    offsets will give us that number anyway, so wait to create the edges.
+   *  - The number of edges is a guess that doesn't include deduplication of new edges with
+   *    existing edges. If those are found, the mesh will be resized later.
+   *  - Don't create attributes to facilite implicit sharing of the positions array. */
+  Mesh *mesh = create_mesh_no_attributes(src_mesh,
+                                         src_mesh.verts_num,
+                                         src_edges.size() + tri_edges_range.size(),
+                                         tris_range.size() + unselected.size(),
+                                         0);
+
+  /* Find the face corner ranges using the offsets array from the new mesh. That gives us the
+   * final number of face corners. */
+  const OffsetIndices faces = calc_faces(src_faces, unselected, mesh->face_offsets_for_write());
+  mesh->corners_num = faces.total_size();
+  const OffsetIndices faces_unselected = faces.slice(unselected_range);
+
+  bke::MutableAttributeAccessor attributes = mesh->attributes_for_write();
+  attributes.add<int2>(".edge_verts", bke::AttrDomain::Edge, bke::AttributeInitConstruct());
+  attributes.add<int>(".corner_vert", bke::AttrDomain::Corner, bke::AttributeInitConstruct());
+  attributes.add<int>(".corner_edge", bke::AttrDomain::Corner, bke::AttributeInitConstruct());
+
+  MutableSpan<int2> edges_with_duplicates = mesh->edges_for_write();
+  MutableSpan<int> corner_verts = mesh->corner_verts_for_write();
+  MutableSpan<int> corner_edges = mesh->corner_edges_for_write();
+
+  bke::attribute_math::gather(src_corner_verts, corner_map, corner_verts.slice(tri_corners_range));
+
+  if (!ngons.is_empty()) {
+    ngon::calc_edges(src_faces,
+                     src_corner_verts,
+                     src_corner_edges,
+                     ngons,
+                     tris_by_ngon,
+                     edges_by_ngon,
+                     faces.slice(ngon_tris_range),
+                     ngon_edges_range.start(),
+                     corner_map.as_mutable_span().slice(ngon_corners_range),
+                     edges_with_duplicates.slice(ngon_edges_range),
+                     corner_edges);
+  }
+
+  if (!quads.is_empty()) {
+    quad::calc_edges(src_corner_edges,
+                     corner_map.as_mutable_span().slice(quad_corners_range),
+                     corner_verts,
+                     quad_edges_range.start(),
+                     edges_with_duplicates.slice(quad_edges_range),
+                     corner_edges.slice(quad_corners_range));
+  }
+
+  {
+    Array<int> vert_to_edge_offsets;
+    Array<int> vert_to_edge_indices;
+    const GroupedSpan<int> src_vert_to_edge_map = bke::mesh::build_vert_to_edge_map(
+        src_edges, src_mesh.verts_num, vert_to_edge_offsets, vert_to_edge_indices);
+
+    mesh->edges_num = deduplication::calc_new_edges(src_mesh,
+                                                    src_edges,
+                                                    src_vert_to_edge_map,
+                                                    tri_edges_range,
+                                                    edges_with_duplicates,
+                                                    corner_edges);
+
+    edges_with_duplicates.take_front(src_edges.size()).copy_from(src_edges);
+  }
+
+  /* Vertex attributes are totally unnaffected and can be shared with implicit sharing.
+   * Use the #CustomData API for better support for vertex groups. */
+  CustomData_merge(&src_mesh.vert_data, &mesh->vert_data, CD_MASK_MESH.vmask, mesh->verts_num);
+
+  for (auto &attribute : bke::retrieve_attributes_for_transfer(
+           src_attributes, attributes, ATTR_DOMAIN_MASK_EDGE, propagation_info, {".edge_verts"}))
+  {
+    attribute.dst.span.slice(src_edges_range).copy_from(attribute.src);
+    GMutableSpan new_data = attribute.dst.span.drop_front(src_edges.size());
+    /* It would be reasonable interpolate data from connected edges within each face.
+     * Currently the data from new edges is just set to the type's default value. */
+    const void *default_value = new_data.type().default_value();
+    new_data.type().fill_construct_n(default_value, new_data.data(), new_data.size());
+    attribute.dst.finish();
+  }
+
+  if (CustomData_has_layer(&src_mesh.edge_data, CD_ORIGINDEX)) {
+    const Span src(
+        static_cast<const int *>(CustomData_get_layer(&src_mesh.edge_data, CD_ORIGINDEX)),
+        src_mesh.edges_num);
+    MutableSpan dst(static_cast<int *>(CustomData_add_layer(
+                        &mesh->edge_data, CD_ORIGINDEX, CD_CONSTRUCT, mesh->edges_num)),
+                    mesh->edges_num);
+    dst.drop_front(src_edges.size()).fill(ORIGINDEX_NONE);
+    array_utils::copy(src, dst.slice(src_edges_range));
+  }
+
+  for (auto &attribute : bke::retrieve_attributes_for_transfer(
+           src_attributes, attributes, ATTR_DOMAIN_MASK_FACE, propagation_info, {}))
+  {
+    bke::attribute_math::gather_to_groups(
+        tris_by_ngon, ngons, attribute.src, attribute.dst.span.slice(ngon_tris_range));
+    quad::copy_face_data_to_tris(attribute.src, quads, attribute.dst.span.slice(quad_tris_range));
+    array_utils::gather(attribute.src, unselected, attribute.dst.span.slice(unselected_range));
+    attribute.dst.finish();
+  }
+
+  if (CustomData_has_layer(&src_mesh.face_data, CD_ORIGINDEX)) {
+    const Span src(
+        static_cast<const int *>(CustomData_get_layer(&src_mesh.face_data, CD_ORIGINDEX)),
+        src_mesh.faces_num);
+    MutableSpan dst(static_cast<int *>(CustomData_add_layer(
+                        &mesh->face_data, CD_ORIGINDEX, CD_CONSTRUCT, mesh->faces_num)),
+                    mesh->faces_num);
+    bke::attribute_math::gather_to_groups(tris_by_ngon, ngons, src, dst.slice(ngon_tris_range));
+    quad::copy_face_data_to_tris(src, quads, dst.slice(quad_tris_range));
+    array_utils::gather(src, unselected, dst.slice(unselected_range));
+  }
+
+  array_utils::gather_group_to_group(
+      src_faces, faces_unselected, unselected, src_corner_verts, corner_verts);
+  array_utils::gather_group_to_group(
+      src_faces, faces_unselected, unselected, src_corner_edges, corner_edges);
+  for (auto &attribute : bke::retrieve_attributes_for_transfer(src_attributes,
+                                                               attributes,
+                                                               ATTR_DOMAIN_MASK_CORNER,
+                                                               propagation_info,
+                                                               {".corner_vert", ".corner_edge"}))
+  {
+    bke::attribute_math::gather_group_to_group(
+        src_faces, faces_unselected, unselected, attribute.src, attribute.dst.span);
+    bke::attribute_math::gather(
+        attribute.src, corner_map.as_span(), attribute.dst.span.slice(tri_corners_range));
+    attribute.dst.finish();
+  }
+
+  mesh->runtime->bounds_cache = src_mesh.runtime->bounds_cache;
+  copy_loose_vert_hint(src_mesh, *mesh);
+  copy_loose_edge_hint(src_mesh, *mesh);
+  if (src_mesh.no_overlapping_topology()) {
+    mesh->tag_overlapping_none();
+  }
+  BLI_assert(BKE_mesh_is_valid(mesh));
+  return mesh;
+}
+
+}  // namespace blender::geometry

source/blender/makesdna/DNA_node_types.h

@@ -2608,19 +2608,6 @@ typedef enum GeometryNodeCurveHandleMode {
   GEO_NODE_CURVE_HANDLE_RIGHT = (1 << 1)
 } GeometryNodeCurveHandleMode;
 
-typedef enum GeometryNodeTriangulateNGons {
-  GEO_NODE_TRIANGULATE_NGON_BEAUTY = 0,
-  GEO_NODE_TRIANGULATE_NGON_EARCLIP = 1,
-} GeometryNodeTriangulateNGons;
-
-typedef enum GeometryNodeTriangulateQuads {
-  GEO_NODE_TRIANGULATE_QUAD_BEAUTY = 0,
-  GEO_NODE_TRIANGULATE_QUAD_FIXED = 1,
-  GEO_NODE_TRIANGULATE_QUAD_ALTERNATE = 2,
-  GEO_NODE_TRIANGULATE_QUAD_SHORTEDGE = 3,
-  GEO_NODE_TRIANGULATE_QUAD_LONGEDGE = 4,
-} GeometryNodeTriangulateQuads;
-
 typedef enum GeometryNodeDistributePointsInVolumeMode {
   GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_RANDOM = 0,
   GEO_NODE_DISTRIBUTE_POINTS_IN_VOLUME_DENSITY_GRID = 1,

source/blender/nodes/geometry/nodes/node_geo_triangulate.cc

@@ -2,16 +2,12 @@
  *
  * SPDX-License-Identifier: GPL-2.0-or-later */
 
-#include "BKE_customdata.hh"
 #include "BKE_mesh.hh"
 
-#include "bmesh.hh"
-#include "bmesh_tools.hh"
-
-#include "DNA_mesh_types.h"
-
 #include "NOD_rna_define.hh"
 
+#include "GEO_mesh_triangulate.hh"
+
 #include "UI_interface.hh"
 #include "UI_resources.hh"
 
@@ -25,7 +21,6 @@ static void node_declare(NodeDeclarationBuilder &b)
 {
   b.add_input<decl::Geometry>("Mesh").supported_type(GeometryComponent::Type::Mesh);
   b.add_input<decl::Bool>("Selection").default_value(true).field_on_all().hide_value();
-  b.add_input<decl::Int>("Minimum Vertices").default_value(4).min(4).max(10000);
   b.add_output<decl::Geometry>("Mesh").propagate_all();
 }
 
@@ -37,69 +32,50 @@ static void node_layout(uiLayout *layout, bContext * /*C*/, PointerRNA *ptr)
 
 static void geo_triangulate_init(bNodeTree * /*tree*/, bNode *node)
 {
-  node->custom1 = GEO_NODE_TRIANGULATE_QUAD_SHORTEDGE;
-  node->custom2 = GEO_NODE_TRIANGULATE_NGON_BEAUTY;
-}
-
-static Mesh *triangulate_mesh_selection(const Mesh &mesh,
-                                        const int quad_method,
-                                        const int ngon_method,
-                                        const IndexMask &selection,
-                                        const int min_vertices)
-{
-  CustomData_MeshMasks cd_mask_extra = {
-      CD_MASK_ORIGINDEX, CD_MASK_ORIGINDEX, 0, CD_MASK_ORIGINDEX};
-  BMeshCreateParams create_params{false};
-  BMeshFromMeshParams from_mesh_params{};
-  from_mesh_params.calc_face_normal = true;
-  from_mesh_params.calc_vert_normal = true;
-  from_mesh_params.cd_mask_extra = cd_mask_extra;
-  BMesh *bm = BKE_mesh_to_bmesh_ex(&mesh, &create_params, &from_mesh_params);
-
-  /* Tag faces to be triangulated from the selection mask. */
-  BM_mesh_elem_table_ensure(bm, BM_FACE);
-  selection.foreach_index([&](const int i_face) {
-    BM_elem_flag_set(BM_face_at_index(bm, i_face), BM_ELEM_TAG, true);
-  });
-
-  BM_mesh_triangulate(bm, quad_method, ngon_method, min_vertices, true, nullptr, nullptr, nullptr);
-  Mesh *result = BKE_mesh_from_bmesh_for_eval_nomain(bm, &cd_mask_extra, &mesh);
-  BM_mesh_free(bm);
-
-  /* Positions are not changed by the triangulation operation, so the bounds are the same. */
-  result->runtime->bounds_cache = mesh.runtime->bounds_cache;
-
-  /* Vertex order is not affected. */
-  geometry::debug_randomize_edge_order(result);
-  geometry::debug_randomize_face_order(result);
-
-  return result;
+  node->custom1 = int(geometry::TriangulateQuadMode::ShortEdge);
+  node->custom2 = int(geometry::TriangulateNGonMode::Beauty);
 }
 
 static void node_geo_exec(GeoNodeExecParams params)
 {
   GeometrySet geometry_set = params.extract_input<GeometrySet>("Mesh");
   Field<bool> selection_field = params.extract_input<Field<bool>>("Selection");
-  const int min_vertices = std::max(params.extract_input<int>("Minimum Vertices"), 4);
+  const AnonymousAttributePropagationInfo &propagation_info = params.get_output_propagation_info(
+      "Mesh");
 
-  GeometryNodeTriangulateQuads quad_method = GeometryNodeTriangulateQuads(params.node().custom1);
-  GeometryNodeTriangulateNGons ngon_method = GeometryNodeTriangulateNGons(params.node().custom2);
+  geometry::TriangulateNGonMode ngon_method = geometry::TriangulateNGonMode(params.node().custom2);
+  geometry::TriangulateQuadMode quad_method = geometry::TriangulateQuadMode(params.node().custom1);
 
   geometry_set.modify_geometry_sets([&](GeometrySet &geometry_set) {
-    if (!geometry_set.has_mesh()) {
+    const Mesh *src_mesh = geometry_set.get_mesh();
+    if (!src_mesh) {
       return;
     }
-    const Mesh &mesh_in = *geometry_set.get_mesh();
 
-    const bke::MeshFieldContext context{mesh_in, AttrDomain::Face};
-    FieldEvaluator evaluator{context, mesh_in.faces_num};
+    const bke::MeshFieldContext context(*src_mesh, AttrDomain::Face);
+    FieldEvaluator evaluator{context, src_mesh->faces_num};
     evaluator.add(selection_field);
     evaluator.evaluate();
     const IndexMask selection = evaluator.get_evaluated_as_mask(0);
+    if (selection.is_empty()) {
+      return;
+    }
 
-    Mesh *mesh_out = triangulate_mesh_selection(
-        mesh_in, quad_method, ngon_method, selection, min_vertices);
-    geometry_set.replace_mesh(mesh_out);
+    std::optional<Mesh *> mesh = geometry::mesh_triangulate(
+        *src_mesh,
+        selection,
+        geometry::TriangulateNGonMode(ngon_method),
+        geometry::TriangulateQuadMode(quad_method),
+        propagation_info);
+    if (!mesh) {
+      return;
+    }
+
+    /* Vertex order is not affected. */
+    geometry::debug_randomize_edge_order(*mesh);
+    geometry::debug_randomize_face_order(*mesh);
+
+    geometry_set.replace_mesh(*mesh);
   });
 
   params.set_output("Mesh", std::move(geometry_set));
@@ -108,27 +84,27 @@ static void node_geo_exec(GeoNodeExecParams params)
 static void node_rna(StructRNA *srna)
 {
   static const EnumPropertyItem rna_node_geometry_triangulate_quad_method_items[] = {
-      {GEO_NODE_TRIANGULATE_QUAD_BEAUTY,
+      {int(geometry::TriangulateQuadMode::Beauty),
        "BEAUTY",
        0,
        "Beauty",
        "Split the quads in nice triangles, slower method"},
-      {GEO_NODE_TRIANGULATE_QUAD_FIXED,
+      {int(geometry::TriangulateQuadMode::Fixed),
        "FIXED",
        0,
        "Fixed",
        "Split the quads on the first and third vertices"},
-      {GEO_NODE_TRIANGULATE_QUAD_ALTERNATE,
+      {int(geometry::TriangulateQuadMode::Alternate),
        "FIXED_ALTERNATE",
        0,
        "Fixed Alternate",
        "Split the quads on the 2nd and 4th vertices"},
-      {GEO_NODE_TRIANGULATE_QUAD_SHORTEDGE,
+      {int(geometry::TriangulateQuadMode::ShortEdge),
        "SHORTEST_DIAGONAL",
        0,
        "Shortest Diagonal",
        "Split the quads along their shortest diagonal"},
-      {GEO_NODE_TRIANGULATE_QUAD_LONGEDGE,
+      {int(geometry::TriangulateQuadMode::LongEdge),
        "LONGEST_DIAGONAL",
        0,
        "Longest Diagonal",
@@ -137,12 +113,12 @@ static void node_rna(StructRNA *srna)
   };
 
   static const EnumPropertyItem rna_node_geometry_triangulate_ngon_method_items[] = {
-      {GEO_NODE_TRIANGULATE_NGON_BEAUTY,
+      {int(geometry::TriangulateNGonMode::Beauty),
        "BEAUTY",
        0,
        "Beauty",
        "Arrange the new triangles evenly (slow)"},
-      {GEO_NODE_TRIANGULATE_NGON_EARCLIP,
+      {int(geometry::TriangulateNGonMode::EarClip),
        "CLIP",
        0,
        "Clip",

tests/python/CMakeLists.txt

@@ -816,6 +816,7 @@ set(geo_node_tests
   mesh
   mesh/extrude
   mesh/split_edges
+  mesh/triangulate
   points
   texture
   utilities