bmo_subdivide.cc

Go to the documentation of this file.

/* SPDX-FileCopyrightText: 2023 Blender Authors
 *
 * SPDX-License-Identifier: GPL-2.0-or-later */

 
#include <array>
 
#include "MEM_guardedalloc.h"
 
#include "BLI_math_geom.h"
#include "BLI_math_vector.h"
#include "BLI_noise.h"
#include "BLI_rand.h"
#include "BLI_stack.h"
#include "BLI_vector.hh"
 
#include "BKE_customdata.hh"
 
#include "bmesh.hh"
#include "intern/bmesh_operators_private.hh"
#include "intern/bmesh_private.hh"
 
using blender::Vector;
 
struct SubDParams {
  int numcuts;
  float smooth;
  int smooth_falloff;
  float fractal;
  float along_normal;
  // int beauty;
  bool use_smooth;
  bool use_smooth_even;
  bool use_sphere;
  bool use_fractal;
  int seed;
  BMOperator *op;
  BMOpSlot *slot_edge_percents;   /* `BMO_slot_get(params->op->slots_in, "edge_percents")`. */
  BMOpSlot *slot_custom_patterns; /* `BMO_slot_get(params->op->slots_in, "custom_patterns")`. */
  float fractal_ofs[3];
 
  /* Runtime storage for shape key. */
  struct {
    int cd_vert_shape_offset;
    int cd_vert_shape_offset_tmp;
    int totlayer;
 
    /* Shape-key holding displaced vertex coordinates for current geometry. */
    int tmpkey;
  } shape_info;
};
 
static void bmo_subd_init_shape_info(BMesh *bm, SubDParams *params)
{
  const int skey = CustomData_number_of_layers(&bm->vdata, CD_SHAPEKEY) - 1;
  params->shape_info.tmpkey = skey;
  params->shape_info.cd_vert_shape_offset = CustomData_get_offset(&bm->vdata, CD_SHAPEKEY);
  params->shape_info.cd_vert_shape_offset_tmp = CustomData_get_n_offset(
      &bm->vdata, CD_SHAPEKEY, skey);
  params->shape_info.totlayer = CustomData_number_of_layers(&bm->vdata, CD_SHAPEKEY);
}
 
using subd_pattern_fill_fp = void (*)(BMesh *bm,
                                      BMFace *face,
                                      BMVert **verts,
                                      const SubDParams *params);
 
/*
 * NOTE: this is a pattern-based edge subdivider.
 * it tries to match a pattern to edge selections on faces,
 * then executes functions to cut them.
 */
struct SubDPattern {
  int seledges[20]; /* selected edges mask, for splitting */
 
  /* verts starts at the first new vert cut, not the first vert in the face */
  subd_pattern_fill_fp connectexec;
  int len; /* total number of verts, before any subdivision */
};
 
/* generic subdivision rules:
 *
 * - two selected edges in a face should make a link
 *   between them.
 *
 * - one edge should do, what? make pretty topology, or just
 *   split the edge only?
 */

#define SUBD_SPLIT 1
 
#define EDGE_PERCENT 2
 
/* I don't think new faces are flagged, currently, but
 * better safe than sorry. */
#define FACE_CUSTOMFILL 4
#define ELE_INNER 8
#define ELE_SPLIT 16

#define FLT_FACE_SPLIT_EPSILON 0.00005f
 
/*
 * NOTE: beauty has been renamed to flag!
 */
 
/* generic subdivision rules:
 *
 * - two selected edges in a face should make a link
 *   between them.
 *
 * - one edge should do, what? make pretty topology, or just
 *   split the edge only?
 */
 
/* connects face with smallest len, which I think should always be correct for
 * edge subdivision */
static BMEdge *connect_smallest_face(BMesh *bm, BMVert *v_a, BMVert *v_b, BMFace **r_f_new)
{
  BMLoop *l_a, *l_b;
  BMFace *f;
 
  /* this isn't the best thing in the world.  it doesn't handle cases where there's
   * multiple faces yet.  that might require a convexity test to figure out which
   * face is "best" and who knows what for non-manifold conditions.
   *
   * NOTE: we allow adjacent here, since there's no chance this happens.
   */
  f = BM_vert_pair_share_face_by_len(v_a, v_b, &l_a, &l_b, true);
 
  if (f) {
    BMFace *f_new;
    BMLoop *l_new;
 
    BLI_assert(!BM_loop_is_adjacent(l_a, l_b));
 
    f_new = BM_face_split(bm, f, l_a, l_b, &l_new, nullptr, false);
 
    if (r_f_new) {
      *r_f_new = f_new;
    }
    return l_new ? l_new->e : nullptr;
  }
 
  return nullptr;
}

static void interp_slerp_co_no_v3(
    const float co_a[3],
    const float no_a[3],
    const float co_b[3],
    const float no_b[3],
    const float no_dir[3], /* caller already knows, avoid normalize */
    float fac,
    float r_co[3])
{
  /* center of the sphere defined by both normals */
  float center[3];
 
  BLI_assert(len_squared_v3v3(no_a, no_b) != 0);
 
  /* calculate sphere 'center' */
  {
    /* use point on plane to */
    float no_mid[3], no_ortho[3];
/* pass this as an arg instead */
#if 0
    float no_dir[3];
#endif
 
    add_v3_v3v3(no_mid, no_a, no_b);
    normalize_v3(no_mid);
 
#if 0
    sub_v3_v3v3(no_dir, co_a, co_b);
    normalize_v3(no_dir);
#endif
 
    /* axis of slerp */
    bool center_ok = false;
    cross_v3_v3v3(no_ortho, no_mid, no_dir);
    if (normalize_v3(no_ortho) != 0.0f) {
      float plane_a[4], plane_b[4], plane_c[4];
      float v_a_no_ortho[3], v_b_no_ortho[3];
 
      /* create planes */
      cross_v3_v3v3(v_a_no_ortho, no_ortho, no_a);
      cross_v3_v3v3(v_b_no_ortho, no_ortho, no_b);
      project_v3_plane(v_a_no_ortho, no_ortho, v_a_no_ortho);
      project_v3_plane(v_b_no_ortho, no_ortho, v_b_no_ortho);
 
      plane_from_point_normal_v3(plane_a, co_a, v_a_no_ortho);
      plane_from_point_normal_v3(plane_b, co_b, v_b_no_ortho);
      plane_from_point_normal_v3(plane_c, co_b, no_ortho);
 
      /* find the sphere center from 3 planes */
      if (isect_plane_plane_plane_v3(plane_a, plane_b, plane_c, center)) {
        center_ok = true;
      }
    }
    if (center_ok == false) {
      mid_v3_v3v3(center, co_a, co_b);
    }
  }
 
  /* calculate the final output 'r_co' */
  {
    float ofs_a[3], ofs_b[3], ofs_slerp[3];
    float dist_a, dist_b;
 
    sub_v3_v3v3(ofs_a, co_a, center);
    sub_v3_v3v3(ofs_b, co_b, center);
 
    dist_a = normalize_v3(ofs_a);
    dist_b = normalize_v3(ofs_b);
 
    if (interp_v3_v3v3_slerp(ofs_slerp, ofs_a, ofs_b, fac)) {
      madd_v3_v3v3fl(r_co, center, ofs_slerp, interpf(dist_b, dist_a, fac));
    }
    else {
      interp_v3_v3v3(r_co, co_a, co_b, fac);
    }
  }
}
 
/* Calculates offset for co, based on fractal, sphere or smooth settings. */
static void alter_co(BMVert *v,
                     BMEdge * /*e_orig*/,
                     const SubDParams *params,
                     const float perc,
                     const BMVert *v_a,
                     const BMVert *v_b)
{
  float *co = static_cast<float *>(
      BM_ELEM_CD_GET_VOID_P(v, params->shape_info.cd_vert_shape_offset_tmp));
  int i;
 
  copy_v3_v3(co, v->co);
 
  if (UNLIKELY(params->use_sphere)) { /* subdivide sphere */
    normalize_v3_length(co, params->smooth);
  }
  else if (params->use_smooth) {
/* calculating twice and blending gives smoother results,
 * removing visible seams. */
#define USE_SPHERE_DUAL_BLEND
 
    const float eps_unit_vec = 1e-5f;
    float smooth;
    float no_dir[3];
 
#ifdef USE_SPHERE_DUAL_BLEND
    float no_reflect[3], co_a[3], co_b[3];
#endif
 
    sub_v3_v3v3(no_dir, v_a->co, v_b->co);
    normalize_v3(no_dir);
 
#ifndef USE_SPHERE_DUAL_BLEND
    if (len_squared_v3v3(v_a->no, v_b->no) < eps_unit_vec) {
      interp_v3_v3v3(co, v_a->co, v_b->co, perc);
    }
    else {
      interp_slerp_co_no_v3(v_a->co, v_a->no, v_b->co, v_b->no, no_dir, perc, co);
    }
#else
    /* sphere-a */
    reflect_v3_v3v3(no_reflect, v_a->no, no_dir);
    if (len_squared_v3v3(v_a->no, no_reflect) < eps_unit_vec) {
      interp_v3_v3v3(co_a, v_a->co, v_b->co, perc);
    }
    else {
      interp_slerp_co_no_v3(v_a->co, v_a->no, v_b->co, no_reflect, no_dir, perc, co_a);
    }
 
    /* sphere-b */
    reflect_v3_v3v3(no_reflect, v_b->no, no_dir);
    if (len_squared_v3v3(v_b->no, no_reflect) < eps_unit_vec) {
      interp_v3_v3v3(co_b, v_a->co, v_b->co, perc);
    }
    else {
      interp_slerp_co_no_v3(v_a->co, no_reflect, v_b->co, v_b->no, no_dir, perc, co_b);
    }
 
    /* blend both spheres */
    interp_v3_v3v3(co, co_a, co_b, perc);
#endif /* USE_SPHERE_DUAL_BLEND */
 
    /* apply falloff */
    if (params->smooth_falloff == SUBD_FALLOFF_LIN) {
      smooth = 1.0f;
    }
    else {
      smooth = fabsf(1.0f - 2.0f * fabsf(0.5f - perc));
      smooth = 1.0f + bmesh_subd_falloff_calc(params->smooth_falloff, smooth);
    }
 
    if (params->use_smooth_even) {
      smooth *= shell_v3v3_mid_normalized_to_dist(v_a->no, v_b->no);
    }
 
    smooth *= params->smooth;
    if (smooth != 1.0f) {
      float co_flat[3];
      interp_v3_v3v3(co_flat, v_a->co, v_b->co, perc);
      interp_v3_v3v3(co, co_flat, co, smooth);
    }
 
#undef USE_SPHERE_DUAL_BLEND
  }
 
  if (params->use_fractal) {
    float normal[3], co2[3], base1[3], base2[3], tvec[3];
    const float len = len_v3v3(v_a->co, v_b->co);
    float fac;
 
    fac = params->fractal * len;
 
    mid_v3_v3v3(normal, v_a->no, v_b->no);
    ortho_basis_v3v3_v3(base1, base2, normal);
 
    add_v3_v3v3(co2, v->co, params->fractal_ofs);
    mul_v3_fl(co2, 10.0f);
 
    tvec[0] = fac *
              (BLI_noise_generic_turbulence(1.0, co2[0], co2[1], co2[2], 15, false, 2) - 0.5f);
    tvec[1] = fac *
              (BLI_noise_generic_turbulence(1.0, co2[1], co2[0], co2[2], 15, false, 2) - 0.5f);
    tvec[2] = fac *
              (BLI_noise_generic_turbulence(1.0, co2[1], co2[2], co2[0], 15, false, 2) - 0.5f);
 
    /* add displacement */
    madd_v3_v3fl(co, normal, tvec[0]);
    madd_v3_v3fl(co, base1, tvec[1] * (1.0f - params->along_normal));
    madd_v3_v3fl(co, base2, tvec[2] * (1.0f - params->along_normal));
  }
 
  /* apply the new difference to the rest of the shape keys,
   * note that this doesn't take rotations into account, we _could_ support
   * this by getting the normals and coords for each shape key and
   * re-calculate the smooth value for each but this is quite involved.
   * for now its ok to simply apply the difference IMHO - campbell */
 
  if (params->shape_info.totlayer > 1) {
    float tvec[3];
 
    sub_v3_v3v3(tvec, v->co, co);
 
    /* skip the last layer since its the temp */
    i = params->shape_info.totlayer - 1;
    co = static_cast<float *>(BM_ELEM_CD_GET_VOID_P(v, params->shape_info.cd_vert_shape_offset));
    while (i--) {
      BLI_assert(co != BM_ELEM_CD_GET_VOID_P(v, params->shape_info.cd_vert_shape_offset_tmp));
      sub_v3_v3(co += 3, tvec);
    }
  }
}
 
/* assumes in the edge is the correct interpolated vertices already */
/* percent defines the interpolation, rad and flag are for special options */
/* results in new vertex with correct coordinate, vertex normal and weight group info */
static BMVert *bm_subdivide_edge_addvert(BMesh *bm,
                                         BMEdge *edge,
                                         BMEdge *e_orig,
                                         const SubDParams *params,
                                         const float factor_edge_split,
                                         const float factor_subd,
                                         BMVert *v_a,
                                         BMVert *v_b,
                                         BMEdge **r_edge)
{
  BMVert *v_new;
 
  v_new = BM_edge_split(bm, edge, edge->v1, r_edge, factor_edge_split);
 
  BMO_vert_flag_enable(bm, v_new, ELE_INNER);
 
  /* offset for smooth or sphere or fractal */
  alter_co(v_new, e_orig, params, factor_subd, v_a, v_b);
 
#if 0  // BMESH_TODO
  /* clip if needed by mirror modifier */
  if (edge->v1->f2) {
    if (edge->v1->f2 & edge->v2->f2 & 1) {
      co[0] = 0.0f;
    }
    if (edge->v1->f2 & edge->v2->f2 & 2) {
      co[1] = 0.0f;
    }
    if (edge->v1->f2 & edge->v2->f2 & 4) {
      co[2] = 0.0f;
    }
  }
#endif
 
  interp_v3_v3v3(v_new->no, v_a->no, v_b->no, factor_subd);
  normalize_v3(v_new->no);
 
  return v_new;
}
 
static BMVert *subdivide_edge_num(BMesh *bm,
                                  BMEdge *edge,
                                  BMEdge *e_orig,
                                  int curpoint,
                                  int totpoint,
                                  const SubDParams *params,
                                  BMVert *v_a,
                                  BMVert *v_b,
                                  BMEdge **r_edge)
{
  BMVert *v_new;
  float factor_edge_split, factor_subd;
 
  if (BMO_edge_flag_test(bm, edge, EDGE_PERCENT) && totpoint == 1) {
    factor_edge_split = BMO_slot_map_float_get(params->slot_edge_percents, edge);
    factor_subd = 0.0f;
  }
  else {
    factor_edge_split = 1.0f / float(totpoint + 1 - curpoint);
    factor_subd = float(curpoint + 1) / float(totpoint + 1);
  }
 
  v_new = bm_subdivide_edge_addvert(
      bm, edge, e_orig, params, factor_edge_split, factor_subd, v_a, v_b, r_edge);
  return v_new;
}
 
static void bm_subdivide_multicut(
    BMesh *bm, BMEdge *edge, const SubDParams *params, BMVert *v_a, BMVert *v_b)
{
  BMEdge *eed = edge, *e_new, e_tmp = *edge;
  BMVert *v, v1_tmp = *edge->v1, v2_tmp = *edge->v2, *v1 = edge->v1, *v2 = edge->v2;
  int i, numcuts = params->numcuts;
 
  e_tmp.v1 = &v1_tmp;
  e_tmp.v2 = &v2_tmp;
 
  for (i = 0; i < numcuts; i++) {
    v = subdivide_edge_num(bm, eed, &e_tmp, i, params->numcuts, params, v_a, v_b, &e_new);
 
    BMO_vert_flag_enable(bm, v, SUBD_SPLIT | ELE_SPLIT);
    BMO_edge_flag_enable(bm, eed, SUBD_SPLIT | ELE_SPLIT);
    BMO_edge_flag_enable(bm, e_new, SUBD_SPLIT | ELE_SPLIT);
 
    BM_CHECK_ELEMENT(v);
    if (v->e) {
      BM_CHECK_ELEMENT(v->e);
    }
    if (v->e && v->e->l) {
      BM_CHECK_ELEMENT(v->e->l->f);
    }
  }
 
  alter_co(v1, &e_tmp, params, 0, &v1_tmp, &v2_tmp);
  alter_co(v2, &e_tmp, params, 1.0, &v1_tmp, &v2_tmp);
}
 
/* NOTE: the patterns are rotated as necessary to
 * match the input geometry.  they're based on the
 * pre-split state of the face */

static void quad_1edge_split(BMesh *bm,
                             BMFace * /*face*/,
                             BMVert **verts,
                             const SubDParams *params)
{
  BMFace *f_new;
  int i, add, numcuts = params->numcuts;
 
  /* if it's odd, the middle face is a quad, otherwise it's a triangle */
  if ((numcuts % 2) == 0) {
    add = 2;
    for (i = 0; i < numcuts; i++) {
      if (i == numcuts / 2) {
        add -= 1;
      }
      connect_smallest_face(bm, verts[i], verts[numcuts + add], &f_new);
    }
  }
  else {
    add = 2;
    for (i = 0; i < numcuts; i++) {
      connect_smallest_face(bm, verts[i], verts[numcuts + add], &f_new);
      if (i == numcuts / 2) {
        add -= 1;
        connect_smallest_face(bm, verts[i], verts[numcuts + add], &f_new);
      }
    }
  }
}
 
static const SubDPattern quad_1edge = {
    {1, 0, 0, 0},
    quad_1edge_split,
    4,
};

static void quad_2edge_split_path(BMesh *bm,
                                  BMFace * /*face*/,
                                  BMVert **verts,
                                  const SubDParams *params)
{
  BMFace *f_new;
  int i, numcuts = params->numcuts;
 
  for (i = 0; i < numcuts; i++) {
    connect_smallest_face(bm, verts[i], verts[numcuts + (numcuts - i)], &f_new);
  }
  connect_smallest_face(bm, verts[numcuts * 2 + 3], verts[numcuts * 2 + 1], &f_new);
}
 
static const SubDPattern quad_2edge_path = {
    {1, 1, 0, 0},
    quad_2edge_split_path,
    4,
};

static void quad_2edge_split_innervert(BMesh *bm,
                                       BMFace * /*face*/,
                                       BMVert **verts,
                                       const SubDParams *params)
{
  BMFace *f_new;
  BMVert *v, *v_last;
  BMEdge *e, *e_new, e_tmp;
  int i, numcuts = params->numcuts;
 
  v_last = verts[numcuts];
 
  for (i = numcuts - 1; i >= 0; i--) {
    e = connect_smallest_face(bm, verts[i], verts[numcuts + (numcuts - i)], &f_new);
 
    e_tmp = *e;
    v = bm_subdivide_edge_addvert(bm, e, &e_tmp, params, 0.5f, 0.5f, e->v1, e->v2, &e_new);
 
    if (i != numcuts - 1) {
      connect_smallest_face(bm, v_last, v, &f_new);
    }
 
    v_last = v;
  }
 
  connect_smallest_face(bm, v_last, verts[numcuts * 2 + 2], &f_new);
}
 
static const SubDPattern quad_2edge_innervert = {
    {1, 1, 0, 0},
    quad_2edge_split_innervert,
    4,
};

static void quad_2edge_split_fan(BMesh *bm,
                                 BMFace * /*face*/,
                                 BMVert **verts,
                                 const SubDParams *params)
{
  BMFace *f_new;
  // BMVert *v;                 /* UNUSED */
  // BMVert *v_last = verts[2]; /* UNUSED */
  // BMEdge *e, *e_new;         /* UNUSED */
  int i, numcuts = params->numcuts;
 
  for (i = 0; i < numcuts; i++) {
    connect_smallest_face(bm, verts[i], verts[numcuts * 2 + 2], &f_new);
    connect_smallest_face(bm, verts[numcuts + (numcuts - i)], verts[numcuts * 2 + 2], &f_new);
  }
}
 
static const SubDPattern quad_2edge_fan = {
    {1, 1, 0, 0},
    quad_2edge_split_fan,
    4,
};

static void quad_3edge_split(BMesh *bm,
                             BMFace * /*face*/,
                             BMVert **verts,
                             const SubDParams *params)
{
  BMFace *f_new;
  int i, add = 0, numcuts = params->numcuts;
 
  for (i = 0; i < numcuts; i++) {
    if (i == numcuts / 2) {
      if (numcuts % 2 != 0) {
        connect_smallest_face(bm, verts[numcuts - i - 1 + add], verts[i + numcuts + 1], &f_new);
      }
      add = numcuts * 2 + 2;
    }
    connect_smallest_face(bm, verts[numcuts - i - 1 + add], verts[i + numcuts + 1], &f_new);
  }
 
  for (i = 0; i < numcuts / 2 + 1; i++) {
    connect_smallest_face(bm, verts[i], verts[(numcuts - i) + numcuts * 2 + 1], &f_new);
  }
}
 
static const SubDPattern quad_3edge = {
    {1, 1, 1, 0},
    quad_3edge_split,
    4,
};

static void quad_4edge_subdivide(BMesh *bm,
                                 BMFace * /*face*/,
                                 BMVert **verts,
                                 const SubDParams *params)
{
  BMFace *f_new;
  BMVert *v, *v1, *v2;
  BMEdge *e, *e_new, e_tmp;
  BMVert **lines;
  int numcuts = params->numcuts;
  int i, j, a, b, s = numcuts + 2 /* , totv = numcuts * 4 + 4 */;
 
  lines = MEM_calloc_arrayN<BMVert *>(size_t(numcuts + 2) * size_t(numcuts + 2), "q_4edge_split");
  /* build a 2-dimensional array of verts,
   * containing every vert (and all new ones)
   * in the face */
 
  /* first line */
  for (i = 0; i < numcuts + 2; i++) {
    lines[i] = verts[numcuts * 3 + 2 + (numcuts - i + 1)];
  }
 
  /* last line */
  for (i = 0; i < numcuts + 2; i++) {
    lines[(s - 1) * s + i] = verts[numcuts + i];
  }
 
  /* first and last members of middle lines */
  for (i = 0; i < numcuts; i++) {
    a = i;
    b = numcuts + 1 + numcuts + 1 + (numcuts - i - 1);
 
    e = connect_smallest_face(bm, verts[a], verts[b], &f_new);
    if (!e) {
      continue;
    }
 
    BMO_edge_flag_enable(bm, e, ELE_INNER);
    BMO_face_flag_enable(bm, f_new, ELE_INNER);
 
    v1 = lines[(i + 1) * s] = verts[a];
    v2 = lines[(i + 1) * s + s - 1] = verts[b];
 
    e_tmp = *e;
    for (a = 0; a < numcuts; a++) {
      v = subdivide_edge_num(bm, e, &e_tmp, a, numcuts, params, v1, v2, &e_new);
 
      BMESH_ASSERT(v != nullptr);
 
      BMO_edge_flag_enable(bm, e_new, ELE_INNER);
      lines[(i + 1) * s + a + 1] = v;
    }
  }
 
  for (i = 1; i < numcuts + 2; i++) {
    for (j = 1; j <= numcuts; j++) {
      a = i * s + j;
      b = (i - 1) * s + j;
      e = connect_smallest_face(bm, lines[a], lines[b], &f_new);
      if (!e) {
        continue;
      }
 
      BMO_edge_flag_enable(bm, e, ELE_INNER);
      BMO_face_flag_enable(bm, f_new, ELE_INNER);
    }
  }
 
  MEM_freeN(lines);
}

static void tri_1edge_split(BMesh *bm, BMFace * /*face*/, BMVert **verts, const SubDParams *params)
{
  BMFace *f_new;
  int i, numcuts = params->numcuts;
 
  for (i = 0; i < numcuts; i++) {
    connect_smallest_face(bm, verts[i], verts[numcuts + 1], &f_new);
  }
}
 
static const SubDPattern tri_1edge = {
    {1, 0, 0},
    tri_1edge_split,
    3,
};

static void tri_3edge_subdivide(BMesh *bm,
                                BMFace * /*face*/,
                                BMVert **verts,
                                const SubDParams *params)
{
  BMFace *f_new;
  BMEdge *e, *e_new, e_tmp;
  BMVert ***lines, *v, v1_tmp, v2_tmp;
  void *stackarr[1];
  int i, j, a, b, numcuts = params->numcuts;
 
  /* number of verts in each lin */
  lines = static_cast<BMVert ***>(
      MEM_callocN(sizeof(void *) * (numcuts + 2), "triangle vert table"));
 
  lines[0] = (BMVert **)stackarr;
  lines[0][0] = verts[numcuts * 2 + 1];
 
  lines[numcuts + 1] = MEM_calloc_arrayN<BMVert *>(numcuts + 2, "triangle vert table 2");
  for (i = 0; i < numcuts; i++) {
    lines[numcuts + 1][i + 1] = verts[i];
  }
  lines[numcuts + 1][0] = verts[numcuts * 3 + 2];
  lines[numcuts + 1][numcuts + 1] = verts[numcuts];
 
  for (i = 0; i < numcuts; i++) {
    lines[i + 1] = MEM_calloc_arrayN<BMVert *>(2 + i, "triangle vert table row");
    a = numcuts * 2 + 2 + i;
    b = numcuts + numcuts - i;
    e = connect_smallest_face(bm, verts[a], verts[b], &f_new);
    if (!e) {
      goto cleanup;
    }
 
    BMO_edge_flag_enable(bm, e, ELE_INNER);
    BMO_face_flag_enable(bm, f_new, ELE_INNER);
 
    lines[i + 1][0] = verts[a];
    lines[i + 1][i + 1] = verts[b];
 
    e_tmp = *e;
    v1_tmp = *verts[a];
    v2_tmp = *verts[b];
    e_tmp.v1 = &v1_tmp;
    e_tmp.v2 = &v2_tmp;
    for (j = 0; j < i; j++) {
      v = subdivide_edge_num(bm, e, &e_tmp, j, i, params, verts[a], verts[b], &e_new);
      lines[i + 1][j + 1] = v;
 
      BMO_edge_flag_enable(bm, e_new, ELE_INNER);
    }
  }

  for (i = 1; i <= numcuts; i++) {
    for (j = 0; j < i; j++) {
      e = connect_smallest_face(bm, lines[i][j], lines[i + 1][j + 1], &f_new);
 
      BMO_edge_flag_enable(bm, e, ELE_INNER);
      BMO_face_flag_enable(bm, f_new, ELE_INNER);
 
      e = connect_smallest_face(bm, lines[i][j + 1], lines[i + 1][j + 1], &f_new);
 
      BMO_edge_flag_enable(bm, e, ELE_INNER);
      BMO_face_flag_enable(bm, f_new, ELE_INNER);
    }
  }
 
cleanup:
  for (i = 1; i < numcuts + 2; i++) {
    if (lines[i]) {
      MEM_freeN(lines[i]);
    }
  }
 
  MEM_freeN(lines);
}
 
static const SubDPattern tri_3edge = {
    {1, 1, 1},
    tri_3edge_subdivide,
    3,
};
 
static const SubDPattern quad_4edge = {
    {1, 1, 1, 1},
    quad_4edge_subdivide,
    4,
};
 
static const SubDPattern *patterns[] = {
    nullptr, /* quad single edge pattern is inserted here */
    nullptr, /* quad corner vert pattern is inserted here */
    nullptr, /* tri single edge pattern is inserted here */
    nullptr,
    &quad_3edge,
    nullptr,
};
 
#define PATTERNS_TOT ARRAY_SIZE(patterns)
 
struct SubDFaceData {
  BMVert *start;
  const SubDPattern *pat;
  int totedgesel; /* only used if pat was nullptr, e.g. no pattern was found */
  BMFace *face;
};
 
void bmo_subdivide_edges_exec(BMesh *bm, BMOperator *op)
{
  BMOpSlot *einput;
  const SubDPattern *pat;
  SubDParams params;
  BLI_Stack *facedata;
  BMIter viter, fiter, liter;
  BMVert *v;
  BMEdge *edge;
  BMLoop *l_new, *l;
  BMFace *face;
  float smooth, fractal, along_normal;
  bool use_sphere, use_single_edge, use_grid_fill, use_only_quads;
  int cornertype, seed, i, j, a, b, numcuts, totesel, smooth_falloff;
 
  BMO_slot_buffer_flag_enable(bm, op->slots_in, "edges", BM_EDGE, SUBD_SPLIT);
 
  numcuts = BMO_slot_int_get(op->slots_in, "cuts");
  seed = BMO_slot_int_get(op->slots_in, "seed");
  smooth = BMO_slot_float_get(op->slots_in, "smooth");
  smooth_falloff = BMO_slot_int_get(op->slots_in, "smooth_falloff");
  fractal = BMO_slot_float_get(op->slots_in, "fractal");
  along_normal = BMO_slot_float_get(op->slots_in, "along_normal");
  cornertype = BMO_slot_int_get(op->slots_in, "quad_corner_type");
 
  use_single_edge = BMO_slot_bool_get(op->slots_in, "use_single_edge");
  use_grid_fill = BMO_slot_bool_get(op->slots_in, "use_grid_fill");
  use_only_quads = BMO_slot_bool_get(op->slots_in, "use_only_quads");
  use_sphere = BMO_slot_bool_get(op->slots_in, "use_sphere");
 
  patterns[1] = nullptr;
  /* straight cut is patterns[1] == nullptr */
  switch (cornertype) {
    case SUBD_CORNER_PATH:
      patterns[1] = &quad_2edge_path;
      break;
    case SUBD_CORNER_INNERVERT:
      patterns[1] = &quad_2edge_innervert;
      break;
    case SUBD_CORNER_FAN:
      patterns[1] = &quad_2edge_fan;
      break;
  }
 
  if (use_single_edge) {
    patterns[0] = &quad_1edge;
    patterns[2] = &tri_1edge;
  }
  else {
    patterns[0] = nullptr;
    patterns[2] = nullptr;
  }
 
  if (use_grid_fill) {
    patterns[3] = &quad_4edge;
    patterns[5] = &tri_3edge;
  }
  else {
    patterns[3] = nullptr;
    patterns[5] = nullptr;
  }
 
  /* Add a temporary shape-key layer to store displacements on current geometry. */
  BM_data_layer_add(bm, &bm->vdata, CD_SHAPEKEY);
 
  bmo_subd_init_shape_info(bm, &params);
 
  BM_ITER_MESH (v, &viter, bm, BM_VERTS_OF_MESH) {
    float *co = static_cast<float *>(
        BM_ELEM_CD_GET_VOID_P(v, params.shape_info.cd_vert_shape_offset_tmp));
    copy_v3_v3(co, v->co);
  }
 
  /* first go through and tag edges */
  BMO_slot_buffer_from_enabled_flag(bm, op, op->slots_in, "edges", BM_EDGE, SUBD_SPLIT);
 
  params.numcuts = numcuts;
  params.op = op;
  params.slot_edge_percents = BMO_slot_get(op->slots_in, "edge_percents");
  params.slot_custom_patterns = BMO_slot_get(op->slots_in, "custom_patterns");
  params.smooth = smooth;
  params.smooth_falloff = smooth_falloff;
  params.seed = seed;
  params.fractal = fractal;
  params.along_normal = along_normal;
  params.use_smooth = (smooth != 0.0f);
  params.use_smooth_even = BMO_slot_bool_get(op->slots_in, "use_smooth_even");
  params.use_fractal = (fractal != 0.0f);
  params.use_sphere = use_sphere;
 
  if (params.use_fractal) {
    RNG *rng = BLI_rng_new_srandom(seed);
 
    params.fractal_ofs[0] = BLI_rng_get_float(rng) * 200.0f;
    params.fractal_ofs[1] = BLI_rng_get_float(rng) * 200.0f;
    params.fractal_ofs[2] = BLI_rng_get_float(rng) * 200.0f;
 
    BLI_rng_free(rng);
  }
 
  BMO_slot_map_to_flag(bm, op->slots_in, "custom_patterns", BM_FACE, FACE_CUSTOMFILL);
 
  BMO_slot_map_to_flag(bm, op->slots_in, "edge_percents", BM_EDGE, EDGE_PERCENT);
 
  facedata = BLI_stack_new(sizeof(SubDFaceData), __func__);
 
  Vector<BMVert *, BM_DEFAULT_ITER_STACK_SIZE> verts;
  Vector<BMEdge *, BM_DEFAULT_ITER_STACK_SIZE> edges;
 
  BM_ITER_MESH (face, &fiter, bm, BM_FACES_OF_MESH) {
    BMEdge *e1 = nullptr, *e2 = nullptr;
    float vec1[3], vec2[3];
    bool matched = false;
 
    /* skip non-quads if requested */
    if (use_only_quads && face->len != 4) {
      continue;
    }
 
    /* figure out which pattern to use */
    verts.reinitialize(face->len);
    edges.reinitialize(face->len);
 
    totesel = 0;
    BM_ITER_ELEM_INDEX (l_new, &liter, face, BM_LOOPS_OF_FACE, i) {
      edges[i] = l_new->e;
      verts[i] = l_new->v;
 
      if (BMO_edge_flag_test(bm, edges[i], SUBD_SPLIT)) {
        if (!e1) {
          e1 = edges[i];
        }
        else {
          e2 = edges[i];
        }
 
        totesel++;
      }
    }
 
    /* make sure the two edges have a valid angle to each other */
    if (totesel == 2 && BM_edge_share_vert_check(e1, e2)) {
      sub_v3_v3v3(vec1, e1->v2->co, e1->v1->co);
      sub_v3_v3v3(vec2, e2->v2->co, e2->v1->co);
      normalize_v3(vec1);
      normalize_v3(vec2);
 
      if (fabsf(dot_v3v3(vec1, vec2)) > 1.0f - FLT_FACE_SPLIT_EPSILON) {
        totesel = 0;
      }
    }
 
    if (BMO_face_flag_test(bm, face, FACE_CUSTOMFILL)) {
      pat = static_cast<const SubDPattern *>(
          *BMO_slot_map_data_get(params.slot_custom_patterns, face));
      for (i = 0; i < pat->len; i++) {
        matched = true;
        for (j = 0; j < pat->len; j++) {
          a = (j + i) % pat->len;
          if (!!BMO_edge_flag_test(bm, edges[a], SUBD_SPLIT) != (!!pat->seledges[j])) {
            matched = false;
            break;
          }
        }
        if (matched) {
          SubDFaceData *fd;
 
          fd = static_cast<SubDFaceData *>(BLI_stack_push_r(facedata));
          fd->pat = pat;
          fd->start = verts[i];
          fd->face = face;
          fd->totedgesel = totesel;
          BMO_face_flag_enable(bm, face, SUBD_SPLIT);
          break;
        }
      }
 
      /* Obviously don't test for other patterns matching. */
      continue;
    }
 
    for (i = 0; i < PATTERNS_TOT; i++) {
      pat = patterns[i];
      if (!pat) {
        continue;
      }
 
      if (pat->len == face->len) {
        for (a = 0; a < pat->len; a++) {
          matched = true;
          for (b = 0; b < pat->len; b++) {
            j = (b + a) % pat->len;
            if (!!BMO_edge_flag_test(bm, edges[j], SUBD_SPLIT) != (!!pat->seledges[b])) {
              matched = false;
              break;
            }
          }
          if (matched) {
            break;
          }
        }
        if (matched) {
          SubDFaceData *fd;
 
          BMO_face_flag_enable(bm, face, SUBD_SPLIT);
 
          fd = static_cast<SubDFaceData *>(BLI_stack_push_r(facedata));
          fd->pat = pat;
          fd->start = verts[a];
          fd->face = face;
          fd->totedgesel = totesel;
          break;
        }
      }
    }
 
    if (!matched && totesel) {
      SubDFaceData *fd;
 
      BMO_face_flag_enable(bm, face, SUBD_SPLIT);
 
      /* must initialize all members here */
      fd = static_cast<SubDFaceData *>(BLI_stack_push_r(facedata));
      fd->start = nullptr;
      fd->pat = nullptr;
      fd->totedgesel = totesel;
      fd->face = face;
    }
  }
 
  einput = BMO_slot_get(op->slots_in, "edges");
 
  /* go through and split edges */
  for (i = 0; i < einput->len; i++) {
    edge = static_cast<BMEdge *>(einput->data.buf[i]);
    bm_subdivide_multicut(bm, edge, &params, edge->v1, edge->v2);
  }
 
  /* copy original-geometry displacements to current coordinates */
  BM_ITER_MESH (v, &viter, bm, BM_VERTS_OF_MESH) {
    const float *co = static_cast<const float *>(
        BM_ELEM_CD_GET_VOID_P(v, params.shape_info.cd_vert_shape_offset_tmp));
    copy_v3_v3(v->co, co);
  }
 
  using LoopPair = std::array<BMLoop *, 2>;
  Vector<BMLoop *, BM_DEFAULT_ITER_STACK_SIZE> loops;
  Vector<LoopPair, BM_DEFAULT_ITER_STACK_SIZE> loops_split;
  for (; !BLI_stack_is_empty(facedata); BLI_stack_discard(facedata)) {
    SubDFaceData *fd = static_cast<SubDFaceData *>(BLI_stack_peek(facedata));
 
    face = fd->face;
 
    /* figure out which pattern to use */
    pat = fd->pat;
 
    if (!pat && fd->totedgesel == 2) {
      /* ok, no pattern.  we still may be able to do something */
 
      /* for case of two edges, connecting them shouldn't be too hard */
      loops.reinitialize(face->len);
      BM_ITER_ELEM_INDEX (l, &liter, face, BM_LOOPS_OF_FACE, a) {
        loops[a] = l;
      }
 
      const int64_t vlen = loops.size();
 
      /* find the boundary of one of the split edges */
      for (a = 0; a < vlen; a++) {
        if (!BMO_vert_flag_test(bm, loops[a ? (a - 1) : (vlen - 1)]->v, ELE_INNER) &&
            BMO_vert_flag_test(bm, loops[a]->v, ELE_INNER))
        {
          break;
        }
      }
      /* Failure to break means there is an internal error. */
      BLI_assert(a < vlen);
 
      if (BMO_vert_flag_test(bm, loops[(a + numcuts + 1) % vlen]->v, ELE_INNER)) {
        b = (a + numcuts + 1) % vlen;
      }
      else {
        /* find the boundary of the other edge. */
        for (j = 0; j < vlen; j++) {
          b = (j + a + numcuts + 1) % vlen;
          if (!BMO_vert_flag_test(bm, loops[b == 0 ? vlen - 1 : b - 1]->v, ELE_INNER) &&
              BMO_vert_flag_test(bm, loops[b]->v, ELE_INNER))
          {
            break;
          }
        }
      }
 
      b += numcuts - 1;
 
      loops_split.reinitialize(numcuts);
      for (j = 0; j < numcuts; j++) {
        bool ok = true;
 
        /* Check for special case, see: #32500.
         * This edge pair could be used by more than one face,
         * in this case it used to (2.63), split both faces along the same verts
         * while it could be calculated which face should do the split,
         * it's ambiguous, so in this case we're better off to skip them as exceptional cases
         * and not try to be clever guessing which face to cut up.
         *
         * To avoid this case we need to check:
         * Do the verts of each share a face besides the one we are subdividing,
         *  (but not connect to make an edge of that face).
         */
        {
          BMLoop *other_loop;
          BMIter other_fiter;
          BM_ITER_ELEM (other_loop, &other_fiter, loops[a]->v, BM_LOOPS_OF_VERT) {
            if (other_loop->f != face) {
              if (BM_vert_in_face(loops[b]->v, other_loop->f)) {
                /* we assume that these verts are not making an edge in the face */
                BLI_assert(other_loop->prev->v != loops[a]->v);
                BLI_assert(other_loop->next->v != loops[a]->v);
 
                ok = false;
                break;
              }
            }
          }
        }
 
        if (ok == true) {
          loops_split[j][0] = loops[a];
          loops_split[j][1] = loops[b];
        }
        else {
          loops_split[j][0] = nullptr;
          loops_split[j][1] = nullptr;
        }
 
        b = (b - 1) % vlen;
        a = (a + 1) % vlen;
      }
 
      /* Since these are newly created vertices, we don't need to worry about them being legal,
       * ... though there are some cases we _should_ check for
       * - concave corner of an ngon.
       * - 2 edges being used in 2+ ngons.
       */
      //          BM_face_splits_check_legal(bm, face, loops_split, BLI_array_len(loops_split));
 
      for (const LoopPair &loop_split : loops_split) {
        if (loop_split[0]) {
          BMFace *f_new;
          BLI_assert(BM_edge_exists(loop_split[0]->v, loop_split[1]->v) == nullptr);
          f_new = BM_face_split(bm, face, loop_split[0], loop_split[1], &l_new, nullptr, false);
          if (f_new) {
            BMO_edge_flag_enable(bm, l_new->e, ELE_INNER);
          }
        }
      }
 
      continue;
    }
    if (!pat) {
      continue;
    }
 
    a = 0;
    BM_ITER_ELEM_INDEX (l_new, &liter, face, BM_LOOPS_OF_FACE, j) {
      if (l_new->v == fd->start) {
        a = j + 1;
        break;
      }
    }
 
    verts.reinitialize(face->len);
    BM_ITER_ELEM_INDEX (l_new, &liter, face, BM_LOOPS_OF_FACE, j) {
      b = (j - a + face->len) % face->len;
      verts[b] = l_new->v;
    }
 
    BM_CHECK_ELEMENT(face);
    pat->connectexec(bm, face, verts.data(), &params);
  }
 
  /* copy original-geometry displacements to current coordinates */
  BM_ITER_MESH (v, &viter, bm, BM_VERTS_OF_MESH) {
    const float *co = static_cast<const float *>(
        BM_ELEM_CD_GET_VOID_P(v, params.shape_info.cd_vert_shape_offset_tmp));
    copy_v3_v3(v->co, co);
  }
 
  BM_data_layer_free_n(bm, &bm->vdata, CD_SHAPEKEY, params.shape_info.tmpkey);
 
  BLI_stack_free(facedata);
 
  BMO_slot_buffer_from_enabled_flag(
      bm, op, op->slots_out, "geom_inner.out", BM_ALL_NOLOOP, ELE_INNER);
  BMO_slot_buffer_from_enabled_flag(
      bm, op, op->slots_out, "geom_split.out", BM_ALL_NOLOOP, ELE_SPLIT);
 
  BMO_slot_buffer_from_enabled_flag(
      bm, op, op->slots_out, "geom.out", BM_ALL_NOLOOP, ELE_INNER | ELE_SPLIT | SUBD_SPLIT);
}
 
void BM_mesh_esubdivide(BMesh *bm,
                        const char edge_hflag,
                        const float smooth,
                        const short smooth_falloff,
                        const bool use_smooth_even,
                        const float fractal,
                        const float along_normal,
                        const int numcuts,
                        const int seltype,
                        const int cornertype,
                        const short use_single_edge,
                        const short use_grid_fill,
                        const short use_only_quads,
                        const int seed)
{
  BMOperator op;
 
  /* `use_sphere` isn't exposed here since its only used for new primitives. */
  BMO_op_initf(bm,
               &op,
               BMO_FLAG_DEFAULTS,
               "subdivide_edges edges=%he "
               "smooth=%f smooth_falloff=%i use_smooth_even=%b "
               "fractal=%f along_normal=%f "
               "cuts=%i "
               "quad_corner_type=%i "
               "use_single_edge=%b use_grid_fill=%b "
               "use_only_quads=%b "
               "seed=%i",
               edge_hflag,
               smooth,
               smooth_falloff,
               use_smooth_even,
               fractal,
               along_normal,
               numcuts,
               cornertype,
               use_single_edge,
               use_grid_fill,
               use_only_quads,
               seed);
 
  BMO_op_exec(bm, &op);
 
  switch (seltype) {
    case SUBDIV_SELECT_NONE: {
      break;
    }
    case SUBDIV_SELECT_ORIG: {
      /* set the newly created data to be selected */
      if (edge_hflag & BM_ELEM_SELECT) {
        BMO_slot_buffer_hflag_enable(
            bm, op.slots_out, "geom_inner.out", BM_ALL_NOLOOP, BM_ELEM_SELECT, true);
        BM_mesh_select_flush_from_verts(bm, true);
      }
      break;
    }
    case SUBDIV_SELECT_INNER: {
      if (edge_hflag & BM_ELEM_SELECT) {
        BMO_slot_buffer_hflag_enable(
            bm, op.slots_out, "geom_inner.out", BM_EDGE | BM_VERT, BM_ELEM_SELECT, true);
      }
      break;
    }
    case SUBDIV_SELECT_LOOPCUT: {
      if (edge_hflag & BM_ELEM_SELECT) {
        /* deselect input */
        BM_mesh_elem_hflag_disable_all(bm, BM_ALL_NOLOOP, BM_ELEM_SELECT, false);
        BMO_slot_buffer_hflag_enable(
            bm, op.slots_out, "geom_inner.out", BM_EDGE, BM_ELEM_SELECT, true);
      }
      break;
    }
  }
  if (edge_hflag & BM_ELEM_SELECT) {
    /* TODO(@ideasman42): the current behavior for face selection flushing
     * can cause parts of disconnected UV islands to become selected.
     * This is caused by the underlying geometry becoming selected.
     * while this is not a bug in UV selection uses are likely to
     * notice this problem when dealing with the UV selections.
     * This can be observed after subdividing the default "Suzanne" model
     * when the UV island of one of the ears is selected.
     *
     * We may want to change the resulting selection after a subdivision
     * to avoid this problem occurring. */
 
    if (bm->uv_select_sync_valid) {
      const int cd_loop_uv_offset = CustomData_get_offset(&bm->ldata, CD_PROP_FLOAT2);
      BM_mesh_uvselect_flush_post_subdivide(bm, cd_loop_uv_offset);
    }
  }
 
  BMO_op_finish(bm, &op);
}
 
void bmo_bisect_edges_exec(BMesh *bm, BMOperator *op)
{
  BMOIter siter;
  BMEdge *e;
  SubDParams params = {0};
 
  params.numcuts = BMO_slot_int_get(op->slots_in, "cuts");
  params.op = op;
  params.slot_edge_percents = BMO_slot_get(op->slots_in, "edge_percents");
 
  BM_data_layer_add(bm, &bm->vdata, CD_SHAPEKEY);
 
  bmo_subd_init_shape_info(bm, &params);
 
  /* tag edges in map */
  BMO_slot_map_to_flag(bm, op->slots_in, "edge_percents", BM_EDGE, EDGE_PERCENT);
 
  /* go through and split edges */
  BMO_ITER (e, &siter, op->slots_in, "edges", BM_EDGE) {
    bm_subdivide_multicut(bm, e, &params, e->v1, e->v2);
  }
 
  BMO_slot_buffer_from_enabled_flag(
      bm, op, op->slots_out, "geom_split.out", BM_ALL_NOLOOP, ELE_SPLIT);
 
  BM_data_layer_free_n(bm, &bm->vdata, CD_SHAPEKEY, params.shape_info.tmpkey);
}