bmo_subdivide_edgering.c

Go to the documentation of this file.

/*
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
#include "MEM_guardedalloc.h"
 
#include "BLI_alloca.h"
#include "BLI_listbase.h"
#include "BLI_math.h"
#include "BLI_utildefines.h"
#include "BLI_utildefines_stack.h"
 
#include "BKE_curve.h"
 
#include "bmesh.h"
 
#include "intern/bmesh_operators_private.h" /* own include */
 
#define VERT_SHARED (1 << 0)
 
#define EDGE_RING (1 << 0)
#define EDGE_RIM (1 << 1)
#define EDGE_IN_STACK (1 << 2)
 
#define FACE_OUT (1 << 0)
#define FACE_SHARED (1 << 1)
#define FACE_IN_STACK (1 << 2)
 
/* -------------------------------------------------------------------- */
/* Specialized Utility Funcs */
 
#ifndef NDEBUG
static uint bm_verts_tag_count(BMesh *bm)
{
  int count = 0;
  BMIter iter;
  BMVert *v;
  BM_ITER_MESH (v, &iter, bm, BM_VERTS_OF_MESH) {
    if (BM_elem_flag_test(v, BM_ELEM_TAG)) {
      count++;
    }
  }
  return count;
}
#endif
 
static float bezier_handle_calc_length_v3(const float co_a[3],
                                          const float no_a[3],
                                          const float co_b[3],
                                          const float no_b[3])
{
  const float dot = dot_v3v3(no_a, no_b);
  /* gives closest approx at a circle with 2 parallel handles */
  float fac = 1.333333f;
  float len;
  if (dot < 0.0f) {
    /* scale down to 0.666 if we point directly at each other rough but ok */
    /* TODO, current blend from dot may not be optimal but its also a detail */
    const float t = 1.0f + dot;
    fac = (fac * t) + (0.75f * (1.0f - t));
  }
 
#if 0
  len = len_v3v3(co_a, co_b);
#else
  /* 2d length projected on plane of normals */
  {
    float co_a_ofs[3];
    cross_v3_v3v3(co_a_ofs, no_a, no_b);
    if (len_squared_v3(co_a_ofs) > FLT_EPSILON) {
      add_v3_v3(co_a_ofs, co_a);
      closest_to_line_v3(co_a_ofs, co_b, co_a, co_a_ofs);
    }
    else {
      copy_v3_v3(co_a_ofs, co_a);
    }
    len = len_v3v3(co_a_ofs, co_b);
  }
#endif
 
  return (len * 0.5f) * fac;
}
 
static void bm_edgeloop_vert_tag(struct BMEdgeLoopStore *el_store, const bool tag)
{
  LinkData *node = BM_edgeloop_verts_get(el_store)->first;
  do {
    BM_elem_flag_set((BMVert *)node->data, BM_ELEM_TAG, tag);
  } while ((node = node->next));
}
 
static void bmo_edgeloop_vert_tag(BMesh *bm,
                                  struct BMEdgeLoopStore *el_store,
                                  const short oflag,
                                  const bool tag)
{
  LinkData *node = BM_edgeloop_verts_get(el_store)->first;
  do {
    BMO_vert_flag_set(bm, (BMVert *)node->data, oflag, tag);
  } while ((node = node->next));
}
 
static bool bmo_face_is_vert_tag_all(BMesh *bm, BMFace *f, short oflag)
{
  BMLoop *l_iter, *l_first;
  l_iter = l_first = BM_FACE_FIRST_LOOP(f);
  do {
    if (!BMO_vert_flag_test(bm, l_iter->v, oflag)) {
      return false;
    }
  } while ((l_iter = l_iter->next) != l_first);
  return true;
}
 
static bool bm_vert_is_tag_edge_connect(BMesh *bm, BMVert *v)
{
  BMIter eiter;
  BMEdge *e;
 
  BM_ITER_ELEM (e, &eiter, v, BM_EDGES_OF_VERT) {
    if (BMO_edge_flag_test(bm, e, EDGE_RING)) {
      BMVert *v_other = BM_edge_other_vert(e, v);
      if (BM_elem_flag_test(v_other, BM_ELEM_TAG)) {
        return true;
      }
    }
  }
  return false;
}
 
/* for now we need full overlap,
 * supporting partial overlap could be done but gets complicated
 * when trimming endpoints is not enough to ensure consistency.
 */
static bool bm_edgeloop_check_overlap_all(BMesh *bm,
                                          struct BMEdgeLoopStore *el_store_a,
                                          struct BMEdgeLoopStore *el_store_b)
{
  bool has_overlap = true;
  LinkData *node;
 
  ListBase *lb_a = BM_edgeloop_verts_get(el_store_a);
  ListBase *lb_b = BM_edgeloop_verts_get(el_store_b);
 
  bm_edgeloop_vert_tag(el_store_a, false);
  bm_edgeloop_vert_tag(el_store_b, true);
 
  for (node = lb_a->first; node; node = node->next) {
    if (bm_vert_is_tag_edge_connect(bm, node->data) == false) {
      has_overlap = false;
      goto finally;
    }
  }
 
  bm_edgeloop_vert_tag(el_store_a, true);
  bm_edgeloop_vert_tag(el_store_b, false);
 
  for (node = lb_b->first; node; node = node->next) {
    if (bm_vert_is_tag_edge_connect(bm, node->data) == false) {
      has_overlap = false;
      goto finally;
    }
  }
 
finally:
  bm_edgeloop_vert_tag(el_store_a, false);
  bm_edgeloop_vert_tag(el_store_b, false);
  return has_overlap;
}
 
/* -------------------------------------------------------------------- */
/* Edge Loop Pairs */
/* key (ordered loop pointers) */
static GSet *bm_edgering_pair_calc(BMesh *bm, ListBase *eloops_rim)
{
  GSet *eloop_pair_gs = BLI_gset_pair_new(__func__);
  GHash *vert_eloop_gh = BLI_ghash_ptr_new(__func__);
 
  struct BMEdgeLoopStore *el_store;
 
  /* create vert -> eloop map */
  for (el_store = eloops_rim->first; el_store; el_store = BM_EDGELOOP_NEXT(el_store)) {
    LinkData *node = BM_edgeloop_verts_get(el_store)->first;
    do {
      BLI_ghash_insert(vert_eloop_gh, node->data, el_store);
    } while ((node = node->next));
  }
 
  /* collect eloop pairs */
  for (el_store = eloops_rim->first; el_store; el_store = BM_EDGELOOP_NEXT(el_store)) {
    BMIter eiter;
    BMEdge *e;
 
    BMVert *v = ((LinkData *)BM_edgeloop_verts_get(el_store)->first)->data;
 
    BM_ITER_ELEM (e, &eiter, (BMVert *)v, BM_EDGES_OF_VERT) {
      if (BMO_edge_flag_test(bm, e, EDGE_RING)) {
        struct BMEdgeLoopStore *el_store_other;
        BMVert *v_other = BM_edge_other_vert(e, v);
        GHashPair pair_test;
 
        el_store_other = BLI_ghash_lookup(vert_eloop_gh, v_other);
 
        /* in rare cases we cant find a match */
        if (el_store_other) {
          pair_test.first = el_store;
          pair_test.second = el_store_other;
 
          if (pair_test.first > pair_test.second) {
            SWAP(const void *, pair_test.first, pair_test.second);
          }
 
          void **pair_key_p;
          if (!BLI_gset_ensure_p_ex(eloop_pair_gs, &pair_test, &pair_key_p)) {
            *pair_key_p = BLI_ghashutil_pairalloc(pair_test.first, pair_test.second);
          }
        }
      }
    }
  }
 
  BLI_ghash_free(vert_eloop_gh, NULL, NULL);
 
  if (BLI_gset_len(eloop_pair_gs) == 0) {
    BLI_gset_free(eloop_pair_gs, NULL);
    eloop_pair_gs = NULL;
  }
 
  return eloop_pair_gs;
}
 
/* -------------------------------------------------------------------- */
/* Subdivide an edge 'n' times and return an open edgeloop */
 
static void bm_edge_subdiv_as_loop(
    BMesh *bm, ListBase *eloops, BMEdge *e, BMVert *v_a, const int cuts)
{
  struct BMEdgeLoopStore *eloop;
  BMVert **v_arr = BLI_array_alloca(v_arr, cuts + 2);
  BMVert *v_b;
  BLI_assert(BM_vert_in_edge(e, v_a));
 
  v_b = BM_edge_other_vert(e, v_a);
 
  BM_edge_split_n(bm, e, cuts, &v_arr[1]);
  if (v_a == e->v1) {
    v_arr[0] = v_a;
    v_arr[cuts + 1] = v_b;
  }
  else {
    v_arr[0] = v_b;
    v_arr[cuts + 1] = v_a;
  }
 
  eloop = BM_edgeloop_from_verts(v_arr, cuts + 2, false);
 
  if (v_a == e->v1) {
    BM_edgeloop_flip(bm, eloop);
  }
 
  BLI_addtail(eloops, eloop);
}
 
/* -------------------------------------------------------------------- */
/* LoopPair Cache (struct and util funcs) */
 
static void bm_vert_calc_surface_tangent(BMesh *bm, BMVert *v, float r_no[3])
{
  BMIter eiter;
  BMEdge *e;
 
  /* get outer normal, fallback to inner (if this vertex is on a boundary) */
  bool found_outer = false, found_inner = false, found_outer_tag = false;
 
  float no_outer[3] = {0.0f}, no_inner[3] = {0.0f};
 
  /* first find rim edges, typically we will only add 2 normals */
  BM_ITER_ELEM (e, &eiter, v, BM_EDGES_OF_VERT) {
    if (UNLIKELY(BM_edge_is_wire(e))) {
      /* pass - this may confuse things */
    }
    else if (BMO_edge_flag_test(bm, e, EDGE_RIM)) {
      BMIter liter;
      BMLoop *l;
      BM_ITER_ELEM (l, &liter, e, BM_LOOPS_OF_EDGE) {
        /* use unmarked (surrounding) faces to create surface tangent */
        float no[3];
        // BM_face_normal_update(l->f);
        BM_edge_calc_face_tangent(e, l, no);
        if (BMO_face_flag_test(bm, l->f, FACE_SHARED)) {
          add_v3_v3(no_inner, no);
          found_inner = true;
        }
        else {
          add_v3_v3(no_outer, no);
          found_outer = true;
 
          /* other side is used too, blend midway */
          if (BMO_face_flag_test(bm, l->f, FACE_OUT)) {
            found_outer_tag = true;
          }
        }
      }
    }
  }
 
  /* detect if this vertex is in-between 2 loops (when blending multiple),
   * if so - take both inner and outer into account */
 
  if (found_inner && found_outer_tag) {
    /* blend between the 2 */
    negate_v3(no_outer);
    normalize_v3(no_outer);
    normalize_v3(no_inner);
    add_v3_v3v3(r_no, no_outer, no_inner);
    normalize_v3(r_no);
  }
  else if (found_outer) {
    negate_v3(no_outer);
    normalize_v3_v3(r_no, no_outer);
  }
  else {
    /* we always have inner geometry */
    BLI_assert(found_inner == true);
    normalize_v3_v3(r_no, no_inner);
  }
}
 
static void bm_faces_share_tag_flush(BMesh *bm, BMEdge **e_arr, const uint e_arr_len)
{
  uint i;
 
  for (i = 0; i < e_arr_len; i++) {
    BMEdge *e = e_arr[i];
    BMLoop *l_iter, *l_first;
 
    l_iter = l_first = e->l;
    do {
      if (!BMO_face_flag_test(bm, l_iter->f, FACE_SHARED)) {
        if (bmo_face_is_vert_tag_all(bm, l_iter->f, VERT_SHARED)) {
          BMO_face_flag_enable(bm, l_iter->f, FACE_SHARED);
        }
      }
    } while ((l_iter = l_iter->radial_next) != l_first);
  }
}
 
static void bm_faces_share_tag_clear(BMesh *bm, BMEdge **e_arr_iter, const uint e_arr_len_iter)
{
  uint i;
 
  for (i = 0; i < e_arr_len_iter; i++) {
    BMEdge *e = e_arr_iter[i];
    BMLoop *l_iter, *l_first;
 
    l_iter = l_first = e->l;
    do {
      BMO_face_flag_disable(bm, l_iter->f, FACE_SHARED);
    } while ((l_iter = l_iter->radial_next) != l_first);
  }
}
 
typedef struct LoopPairStore {
  /* handle array for splines */
  float (*nors_a)[3];
  float (*nors_b)[3];
 
  /* since we don't have reliable index values into the array,
   * store a map (BMVert -> index) */
  GHash *nors_gh_a;
  GHash *nors_gh_b;
} LoopPairStore;
 
static LoopPairStore *bm_edgering_pair_store_create(BMesh *bm,
                                                    struct BMEdgeLoopStore *el_store_a,
                                                    struct BMEdgeLoopStore *el_store_b,
                                                    const int interp_mode)
{
  LoopPairStore *lpair = MEM_mallocN(sizeof(*lpair), __func__);
 
  if (interp_mode == SUBD_RING_INTERP_SURF) {
    const uint len_a = BM_edgeloop_length_get(el_store_a);
    const uint len_b = BM_edgeloop_length_get(el_store_b);
    const uint e_arr_a_len = len_a - (BM_edgeloop_is_closed(el_store_a) ? 0 : 1);
    const uint e_arr_b_len = len_b - (BM_edgeloop_is_closed(el_store_b) ? 0 : 1);
    BMEdge **e_arr_a = BLI_array_alloca(e_arr_a, e_arr_a_len);
    BMEdge **e_arr_b = BLI_array_alloca(e_arr_b, e_arr_b_len);
    uint i;
 
    struct BMEdgeLoopStore *el_store_pair[2] = {el_store_a, el_store_b};
    uint side_index;
    float(*nors_pair[2])[3];
    GHash *nors_gh_pair[2];
 
    BM_edgeloop_edges_get(el_store_a, e_arr_a);
    BM_edgeloop_edges_get(el_store_b, e_arr_b);
 
    lpair->nors_a = MEM_mallocN(sizeof(*lpair->nors_a) * len_a, __func__);
    lpair->nors_b = MEM_mallocN(sizeof(*lpair->nors_b) * len_b, __func__);
 
    nors_pair[0] = lpair->nors_a;
    nors_pair[1] = lpair->nors_b;
 
    lpair->nors_gh_a = BLI_ghash_ptr_new(__func__);
    lpair->nors_gh_b = BLI_ghash_ptr_new(__func__);
 
    nors_gh_pair[0] = lpair->nors_gh_a;
    nors_gh_pair[1] = lpair->nors_gh_b;
 
    /* now calculate nor */
 
    /* all other verts must _not_ be tagged */
    bmo_edgeloop_vert_tag(bm, el_store_a, VERT_SHARED, true);
    bmo_edgeloop_vert_tag(bm, el_store_b, VERT_SHARED, true);
 
    /* tag all faces that are in-between both loops */
    bm_faces_share_tag_flush(bm, e_arr_a, e_arr_a_len);
    bm_faces_share_tag_flush(bm, e_arr_b, e_arr_b_len);
 
    /* now we have all data we need, calculate vertex spline nor! */
    for (side_index = 0; side_index < 2; side_index++) {
      /* iter vars */
      struct BMEdgeLoopStore *el_store = el_store_pair[side_index];
      ListBase *lb = BM_edgeloop_verts_get(el_store);
      GHash *nors_gh_iter = nors_gh_pair[side_index];
      float(*nor)[3] = nors_pair[side_index];
 
      LinkData *v_iter;
 
      for (v_iter = lb->first, i = 0; v_iter; v_iter = v_iter->next, i++) {
        BMVert *v = v_iter->data;
        bm_vert_calc_surface_tangent(bm, v, nor[i]);
        BLI_ghash_insert(nors_gh_iter, v, POINTER_FROM_UINT(i));
      }
    }
 
    /* cleanup verts share */
    bmo_edgeloop_vert_tag(bm, el_store_a, VERT_SHARED, false);
    bmo_edgeloop_vert_tag(bm, el_store_b, VERT_SHARED, false);
 
    /* cleanup faces share */
    bm_faces_share_tag_clear(bm, e_arr_a, e_arr_a_len);
    bm_faces_share_tag_clear(bm, e_arr_b, e_arr_b_len);
  }
  return lpair;
}
 
static void bm_edgering_pair_store_free(LoopPairStore *lpair, const int interp_mode)
{
  if (interp_mode == SUBD_RING_INTERP_SURF) {
    MEM_freeN(lpair->nors_a);
    MEM_freeN(lpair->nors_b);
 
    BLI_ghash_free(lpair->nors_gh_a, NULL, NULL);
    BLI_ghash_free(lpair->nors_gh_b, NULL, NULL);
  }
  MEM_freeN(lpair);
}
 
/* -------------------------------------------------------------------- */
/* Interpolation Function */
 
static void bm_edgering_pair_interpolate(BMesh *bm,
                                         LoopPairStore *lpair,
                                         struct BMEdgeLoopStore *el_store_a,
                                         struct BMEdgeLoopStore *el_store_b,
                                         ListBase *eloops_ring,
                                         const int interp_mode,
                                         const int cuts,
                                         const float smooth,
                                         const float *falloff_cache)
{
  const int resolu = cuts + 2;
  const int dims = 3;
  bool is_a_no_valid, is_b_no_valid;
  int i;
 
  float el_store_a_co[3], el_store_b_co[3];
  float el_store_a_no[3], el_store_b_no[3];
 
  struct BMEdgeLoopStore *el_store_ring;
 
  float(*coord_array_main)[3] = NULL;
 
  BM_edgeloop_calc_center(bm, el_store_a);
  BM_edgeloop_calc_center(bm, el_store_b);
 
  is_a_no_valid = BM_edgeloop_calc_normal(bm, el_store_a);
  is_b_no_valid = BM_edgeloop_calc_normal(bm, el_store_b);
 
  copy_v3_v3(el_store_a_co, BM_edgeloop_center_get(el_store_a));
  copy_v3_v3(el_store_b_co, BM_edgeloop_center_get(el_store_b));
 
  /* correct normals need to be flipped to face each other
   * we know both normals point in the same direction so one will need flipping */
  {
    float el_dir[3];
    float no[3];
    sub_v3_v3v3(el_dir, el_store_a_co, el_store_b_co);
    normalize_v3_v3(no, el_dir);
 
    if (is_a_no_valid == false) {
      is_a_no_valid = BM_edgeloop_calc_normal_aligned(bm, el_store_a, no);
    }
    if (is_b_no_valid == false) {
      is_b_no_valid = BM_edgeloop_calc_normal_aligned(bm, el_store_b, no);
    }
    (void)is_a_no_valid, (void)is_b_no_valid;
 
    copy_v3_v3(el_store_a_no, BM_edgeloop_normal_get(el_store_a));
    copy_v3_v3(el_store_b_no, BM_edgeloop_normal_get(el_store_b));
 
    if (dot_v3v3(el_store_a_no, el_dir) > 0.0f) {
      negate_v3(el_store_a_no);
    }
    if (dot_v3v3(el_store_b_no, el_dir) < 0.0f) {
      negate_v3(el_store_b_no);
    }
  }
  /* now normals are correct, don't touch! */
 
  /* calculate the center spline, multiple  */
  if ((interp_mode == SUBD_RING_INTERP_PATH) || falloff_cache) {
    float handle_a[3], handle_b[3];
    float handle_len;
 
    handle_len = bezier_handle_calc_length_v3(
                     el_store_a_co, el_store_a_no, el_store_b_co, el_store_b_no) *
                 smooth;
 
    mul_v3_v3fl(handle_a, el_store_a_no, handle_len);
    mul_v3_v3fl(handle_b, el_store_b_no, handle_len);
 
    add_v3_v3(handle_a, el_store_a_co);
    add_v3_v3(handle_b, el_store_b_co);
 
    coord_array_main = MEM_mallocN(dims * (resolu) * sizeof(float), __func__);
 
    for (i = 0; i < dims; i++) {
      BKE_curve_forward_diff_bezier(el_store_a_co[i],
                                    handle_a[i],
                                    handle_b[i],
                                    el_store_b_co[i],
                                    ((float *)coord_array_main) + i,
                                    resolu - 1,
                                    sizeof(float) * dims);
    }
  }
 
  switch (interp_mode) {
    case SUBD_RING_INTERP_LINEAR: {
      if (falloff_cache) {
        float(*coord_array)[3] = MEM_mallocN(dims * (resolu) * sizeof(float), __func__);
        for (i = 0; i < resolu; i++) {
          interp_v3_v3v3(
              coord_array[i], el_store_a_co, el_store_b_co, (float)i / (float)(resolu - 1));
        }
 
        for (el_store_ring = eloops_ring->first; el_store_ring;
             el_store_ring = BM_EDGELOOP_NEXT(el_store_ring)) {
          ListBase *lb_ring = BM_edgeloop_verts_get(el_store_ring);
          LinkData *v_iter;
 
          for (v_iter = lb_ring->first, i = 0; v_iter; v_iter = v_iter->next, i++) {
            if (i > 0 && i < resolu - 1) {
              /* shape */
              if (falloff_cache) {
                interp_v3_v3v3(((BMVert *)v_iter->data)->co,
                               coord_array[i],
                               ((BMVert *)v_iter->data)->co,
                               falloff_cache[i]);
              }
            }
          }
        }
 
        MEM_freeN(coord_array);
      }
 
      break;
    }
    case SUBD_RING_INTERP_PATH: {
      float(*direction_array)[3] = MEM_mallocN(dims * (resolu) * sizeof(float), __func__);
      float(*quat_array)[4] = MEM_mallocN(resolu * sizeof(*quat_array), __func__);
      float(*tri_array)[3][3] = MEM_mallocN(resolu * sizeof(*tri_array), __func__);
      float(*tri_sta)[3], (*tri_end)[3], (*tri_tmp)[3];
 
      /* very similar to make_bevel_list_3D_minimum_twist */
 
      /* calculate normals */
      copy_v3_v3(direction_array[0], el_store_a_no);
      negate_v3_v3(direction_array[resolu - 1], el_store_b_no);
      for (i = 1; i < resolu - 1; i++) {
        bisect_v3_v3v3v3(direction_array[i],
                         coord_array_main[i - 1],
                         coord_array_main[i],
                         coord_array_main[i + 1]);
      }
 
      vec_to_quat(quat_array[0], direction_array[0], 5, 1);
      normalize_qt(quat_array[0]);
 
      for (i = 1; i < resolu; i++) {
        float angle = angle_normalized_v3v3(direction_array[i - 1], direction_array[i]);
        // BLI_assert(angle < DEG2RADF(90.0f));
        if (angle > 0.0f) { /* otherwise we can keep as is */
          float cross_tmp[3];
          float q[4];
          cross_v3_v3v3(cross_tmp, direction_array[i - 1], direction_array[i]);
          axis_angle_to_quat(q, cross_tmp, angle);
          mul_qt_qtqt(quat_array[i], q, quat_array[i - 1]);
          normalize_qt(quat_array[i]);
        }
        else {
          copy_qt_qt(quat_array[i], quat_array[i - 1]);
        }
      }
 
      /* init base tri */
      for (i = 0; i < resolu; i++) {
        int j;
 
        const float shape_size = falloff_cache ? falloff_cache[i] : 1.0f;
 
        tri_tmp = tri_array[i];
 
        /* create the triangle and transform */
        for (j = 0; j < 3; j++) {
          zero_v3(tri_tmp[j]);
          if (j == 1) {
            tri_tmp[j][0] = shape_size;
          }
          else if (j == 2) {
            tri_tmp[j][1] = shape_size;
          }
          mul_qt_v3(quat_array[i], tri_tmp[j]);
          add_v3_v3(tri_tmp[j], coord_array_main[i]);
        }
      }
 
      tri_sta = tri_array[0];
      tri_end = tri_array[resolu - 1];
 
      for (el_store_ring = eloops_ring->first; el_store_ring;
           el_store_ring = BM_EDGELOOP_NEXT(el_store_ring)) {
        ListBase *lb_ring = BM_edgeloop_verts_get(el_store_ring);
        LinkData *v_iter;
 
        BMVert *v_a = ((LinkData *)lb_ring->first)->data;
        BMVert *v_b = ((LinkData *)lb_ring->last)->data;
 
        /* skip first and last */
        for (v_iter = ((LinkData *)lb_ring->first)->next, i = 1; v_iter != lb_ring->last;
             v_iter = v_iter->next, i++) {
          float co_a[3], co_b[3];
 
          tri_tmp = tri_array[i];
 
          transform_point_by_tri_v3(co_a, v_a->co, UNPACK3(tri_tmp), UNPACK3(tri_sta));
          transform_point_by_tri_v3(co_b, v_b->co, UNPACK3(tri_tmp), UNPACK3(tri_end));
 
          interp_v3_v3v3(((BMVert *)v_iter->data)->co, co_a, co_b, (float)i / (float)(resolu - 1));
        }
      }
 
      MEM_freeN(direction_array);
      MEM_freeN(quat_array);
      MEM_freeN(tri_array);
      break;
    }
    case SUBD_RING_INTERP_SURF: {
      float(*coord_array)[3] = MEM_mallocN(dims * (resolu) * sizeof(float), __func__);
 
      /* calculate a bezier handle per edge ring */
      for (el_store_ring = eloops_ring->first; el_store_ring;
           el_store_ring = BM_EDGELOOP_NEXT(el_store_ring)) {
        ListBase *lb_ring = BM_edgeloop_verts_get(el_store_ring);
        LinkData *v_iter;
 
        BMVert *v_a = ((LinkData *)lb_ring->first)->data;
        BMVert *v_b = ((LinkData *)lb_ring->last)->data;
 
        float co_a[3], no_a[3], handle_a[3], co_b[3], no_b[3], handle_b[3];
        float handle_len;
 
        copy_v3_v3(co_a, v_a->co);
        copy_v3_v3(co_b, v_b->co);
 
        /* don't calculate normals here else we get into feedback loop
         * when subdividing 2+ connected edge rings */
#if 0
        bm_vert_calc_surface_tangent(bm, v_a, no_a);
        bm_vert_calc_surface_tangent(bm, v_b, no_b);
#else
        {
          const uint index_a = POINTER_AS_UINT(BLI_ghash_lookup(lpair->nors_gh_a, v_a));
          const uint index_b = POINTER_AS_UINT(BLI_ghash_lookup(lpair->nors_gh_b, v_b));
 
          BLI_assert(BLI_ghash_haskey(lpair->nors_gh_a, v_a));
          BLI_assert(BLI_ghash_haskey(lpair->nors_gh_b, v_b));
 
          copy_v3_v3(no_a, lpair->nors_a[index_a]);
          copy_v3_v3(no_b, lpair->nors_b[index_b]);
        }
#endif
        handle_len = bezier_handle_calc_length_v3(co_a, no_a, co_b, no_b) * smooth;
 
        mul_v3_v3fl(handle_a, no_a, handle_len);
        mul_v3_v3fl(handle_b, no_b, handle_len);
 
        add_v3_v3(handle_a, co_a);
        add_v3_v3(handle_b, co_b);
 
        for (i = 0; i < dims; i++) {
          BKE_curve_forward_diff_bezier(co_a[i],
                                        handle_a[i],
                                        handle_b[i],
                                        co_b[i],
                                        ((float *)coord_array) + i,
                                        resolu - 1,
                                        sizeof(float) * dims);
        }
 
        /* skip first and last */
        for (v_iter = ((LinkData *)lb_ring->first)->next, i = 1; v_iter != lb_ring->last;
             v_iter = v_iter->next, i++) {
          if (i > 0 && i < resolu - 1) {
            copy_v3_v3(((BMVert *)v_iter->data)->co, coord_array[i]);
 
            /* shape */
            if (falloff_cache) {
              interp_v3_v3v3(((BMVert *)v_iter->data)->co,
                             coord_array_main[i],
                             ((BMVert *)v_iter->data)->co,
                             falloff_cache[i]);
            }
          }
        }
      }
 
      MEM_freeN(coord_array);
 
      break;
    }
  }
 
  if (coord_array_main) {
    MEM_freeN(coord_array_main);
  }
}
 
static void bm_face_slice(BMesh *bm, BMLoop *l, const int cuts)
{
  /* TODO, interpolate edge data */
  BMLoop *l_new = l;
  int i;
 
  for (i = 0; i < cuts; i++) {
    /* no chance of double */
    BM_face_split(bm, l_new->f, l_new->prev, l_new->next->next, &l_new, NULL, false);
    if (l_new == NULL) {
      /* This happens when l_new->prev and l_new->next->next are adjacent. Since
       * this sets l_new to NULL, we cannot continue this for-loop. */
      break;
    }
    if (l_new->f->len < l_new->radial_next->f->len) {
      l_new = l_new->radial_next;
    }
    BMO_face_flag_enable(bm, l_new->f, FACE_OUT);
    BMO_face_flag_enable(bm, l_new->radial_next->f, FACE_OUT);
  }
}
 
static bool bm_edgering_pair_order_is_flipped(BMesh *UNUSED(bm),
                                              struct BMEdgeLoopStore *el_store_a,
                                              struct BMEdgeLoopStore *el_store_b)
{
  ListBase *lb_a = BM_edgeloop_verts_get(el_store_a);
  ListBase *lb_b = BM_edgeloop_verts_get(el_store_b);
 
  LinkData *v_iter_a_first = lb_a->first;
  LinkData *v_iter_b_first = lb_b->first;
 
  LinkData *v_iter_a_step = v_iter_a_first;
  LinkData *v_iter_b_step = v_iter_b_first;
 
  /* we _must_ have same starting edge shared */
  BLI_assert(BM_edge_exists(v_iter_a_first->data, v_iter_b_first->data));
 
  /* step around any fan-faces on both sides */
  do {
    v_iter_a_step = v_iter_a_step->next;
  } while (v_iter_a_step && ((BM_edge_exists(v_iter_a_step->data, v_iter_b_first->data)) ||
                             (BM_edge_exists(v_iter_a_step->data, v_iter_b_first->next->data))));
  do {
    v_iter_b_step = v_iter_b_step->next;
  } while (v_iter_b_step && ((BM_edge_exists(v_iter_b_step->data, v_iter_a_first->data)) ||
                             (BM_edge_exists(v_iter_b_step->data, v_iter_a_first->next->data))));
 
  v_iter_a_step = v_iter_a_step ? v_iter_a_step->prev : lb_a->last;
  v_iter_b_step = v_iter_b_step ? v_iter_b_step->prev : lb_b->last;
 
  return !(BM_edge_exists(v_iter_a_step->data, v_iter_b_step->data) ||
           BM_edge_exists(v_iter_a_first->next->data, v_iter_b_step->data) ||
           BM_edge_exists(v_iter_b_first->next->data, v_iter_a_step->data));
}
 
static void bm_edgering_pair_order(BMesh *bm,
                                   struct BMEdgeLoopStore *el_store_a,
                                   struct BMEdgeLoopStore *el_store_b)
{
  ListBase *lb_a = BM_edgeloop_verts_get(el_store_a);
  ListBase *lb_b = BM_edgeloop_verts_get(el_store_b);
 
  LinkData *node;
 
  bm_edgeloop_vert_tag(el_store_a, false);
  bm_edgeloop_vert_tag(el_store_b, true);
 
  /* before going much further, get ourselves in order
   * - align loops (not strictly necessary but handy)
   * - ensure winding is set for both loops */
  if (BM_edgeloop_is_closed(el_store_a) && BM_edgeloop_is_closed(el_store_b)) {
    BMIter eiter;
    BMEdge *e;
    BMVert *v_other;
 
    node = lb_a->first;
 
    BM_ITER_ELEM (e, &eiter, (BMVert *)node->data, BM_EDGES_OF_VERT) {
      if (BMO_edge_flag_test(bm, e, EDGE_RING)) {
        v_other = BM_edge_other_vert(e, (BMVert *)node->data);
        if (BM_elem_flag_test(v_other, BM_ELEM_TAG)) {
          break;
        }
        v_other = NULL;
      }
    }
    BLI_assert(v_other != NULL);
 
    for (node = lb_b->first; node; node = node->next) {
      if (node->data == v_other) {
        break;
      }
    }
    BLI_assert(node != NULL);
 
    BLI_listbase_rotate_first(lb_b, node);
 
    /* now check we are winding the same way */
    if (bm_edgering_pair_order_is_flipped(bm, el_store_a, el_store_b)) {
      BM_edgeloop_flip(bm, el_store_b);
      /* re-ensure the first node */
      BLI_listbase_rotate_first(lb_b, node);
    }
 
    /* sanity checks that we are aligned & winding now */
    BLI_assert(bm_edgering_pair_order_is_flipped(bm, el_store_a, el_store_b) == false);
  }
  else {
    /* If we don't share and edge - flip. */
    BMEdge *e = BM_edge_exists(((LinkData *)lb_a->first)->data, ((LinkData *)lb_b->first)->data);
    if (e == NULL || !BMO_edge_flag_test(bm, e, EDGE_RING)) {
      BM_edgeloop_flip(bm, el_store_b);
    }
  }
 
  /* for cases with multiple loops */
  bm_edgeloop_vert_tag(el_store_b, false);
}
 
static void bm_edgering_pair_subdiv(BMesh *bm,
                                    struct BMEdgeLoopStore *el_store_a,
                                    struct BMEdgeLoopStore *el_store_b,
                                    ListBase *eloops_ring,
                                    const int cuts)
{
  ListBase *lb_a = BM_edgeloop_verts_get(el_store_a);
  // ListBase *lb_b = BM_edgeloop_verts_get(el_store_b);
  const int stack_max = max_ii(BM_edgeloop_length_get(el_store_a),
                               BM_edgeloop_length_get(el_store_b)) *
                        2;
  BMEdge **edges_ring_arr = BLI_array_alloca(edges_ring_arr, stack_max);
  BMFace **faces_ring_arr = BLI_array_alloca(faces_ring_arr, stack_max);
  STACK_DECLARE(edges_ring_arr);
  STACK_DECLARE(faces_ring_arr);
  struct BMEdgeLoopStore *el_store_ring;
  LinkData *node;
  BMEdge *e;
  BMFace *f;
 
  STACK_INIT(edges_ring_arr, stack_max);
  STACK_INIT(faces_ring_arr, stack_max);
 
  bm_edgeloop_vert_tag(el_store_a, false);
  bm_edgeloop_vert_tag(el_store_b, true);
 
  for (node = lb_a->first; node; node = node->next) {
    BMIter eiter;
 
    BM_ITER_ELEM (e, &eiter, (BMVert *)node->data, BM_EDGES_OF_VERT) {
      if (!BMO_edge_flag_test(bm, e, EDGE_IN_STACK)) {
        BMVert *v_other = BM_edge_other_vert(e, (BMVert *)node->data);
        if (BM_elem_flag_test(v_other, BM_ELEM_TAG)) {
          BMIter fiter;
 
          BMO_edge_flag_enable(bm, e, EDGE_IN_STACK);
          STACK_PUSH(edges_ring_arr, e);
 
          /* add faces to the stack */
          BM_ITER_ELEM (f, &fiter, e, BM_FACES_OF_EDGE) {
            if (BMO_face_flag_test(bm, f, FACE_OUT)) {
              if (!BMO_face_flag_test(bm, f, FACE_IN_STACK)) {
                BMO_face_flag_enable(bm, f, FACE_IN_STACK);
                STACK_PUSH(faces_ring_arr, f);
              }
            }
          }
        }
      }
    }
  }
 
  while ((e = STACK_POP(edges_ring_arr))) {
    /* found opposite edge */
    BMVert *v_other;
 
    BMO_edge_flag_disable(bm, e, EDGE_IN_STACK);
 
    /* unrelated to subdiv, but if we _don't_ clear flag, multiple rings fail */
    BMO_edge_flag_disable(bm, e, EDGE_RING);
 
    v_other = BM_elem_flag_test(e->v1, BM_ELEM_TAG) ? e->v1 : e->v2;
    bm_edge_subdiv_as_loop(bm, eloops_ring, e, v_other, cuts);
  }
 
  while ((f = STACK_POP(faces_ring_arr))) {
    BMLoop *l_iter, *l_first;
 
    BMO_face_flag_disable(bm, f, FACE_IN_STACK);
 
    /* Check each edge of the face */
    l_iter = l_first = BM_FACE_FIRST_LOOP(f);
    do {
      if (BMO_edge_flag_test(bm, l_iter->e, EDGE_RIM)) {
        bm_face_slice(bm, l_iter, cuts);
        break;
      }
    } while ((l_iter = l_iter->next) != l_first);
  }
 
  /* clear tags so subdiv verts don't get tagged too  */
  for (el_store_ring = eloops_ring->first; el_store_ring;
       el_store_ring = BM_EDGELOOP_NEXT(el_store_ring)) {
    bm_edgeloop_vert_tag(el_store_ring, false);
  }
 
  /* cleanup after */
  bm_edgeloop_vert_tag(el_store_b, false);
}
 
static void bm_edgering_pair_ringsubd(BMesh *bm,
                                      LoopPairStore *lpair,
                                      struct BMEdgeLoopStore *el_store_a,
                                      struct BMEdgeLoopStore *el_store_b,
                                      const int interp_mode,
                                      const int cuts,
                                      const float smooth,
                                      const float *falloff_cache)
{
  ListBase eloops_ring = {NULL};
  bm_edgering_pair_order(bm, el_store_a, el_store_b);
  bm_edgering_pair_subdiv(bm, el_store_a, el_store_b, &eloops_ring, cuts);
  bm_edgering_pair_interpolate(
      bm, lpair, el_store_a, el_store_b, &eloops_ring, interp_mode, cuts, smooth, falloff_cache);
  BM_mesh_edgeloops_free(&eloops_ring);
}
 
static bool bm_edge_rim_test_cb(BMEdge *e, void *bm_v)
{
  BMesh *bm = bm_v;
  return BMO_edge_flag_test_bool(bm, e, EDGE_RIM);
}
 
/* keep this operator fast, its used in a modifier */
void bmo_subdivide_edgering_exec(BMesh *bm, BMOperator *op)
{
  ListBase eloops_rim = {NULL};
  BMOIter siter;
  BMEdge *e;
  int count;
  bool changed = false;
 
  const int cuts = BMO_slot_int_get(op->slots_in, "cuts");
  const int interp_mode = BMO_slot_int_get(op->slots_in, "interp_mode");
  const float smooth = BMO_slot_float_get(op->slots_in, "smooth");
  const int resolu = cuts + 2;
 
  /* optional 'shape' */
  const int profile_shape = BMO_slot_int_get(op->slots_in, "profile_shape");
  const float profile_shape_factor = BMO_slot_float_get(op->slots_in, "profile_shape_factor");
  float *falloff_cache = (profile_shape_factor != 0.0f) ?
                             BLI_array_alloca(falloff_cache, cuts + 2) :
                             NULL;
 
  BMO_slot_buffer_flag_enable(bm, op->slots_in, "edges", BM_EDGE, EDGE_RING);
 
  BM_mesh_elem_hflag_disable_all(bm, BM_VERT, BM_ELEM_TAG, false);
 
  /* -------------------------------------------------------------------- */
  /* flag outer edges (loops defined as edges on the bounds of the edge ring) */
 
  BMO_ITER (e, &siter, op->slots_in, "edges", BM_EDGE) {
    BMIter fiter;
    BMFace *f;
 
    BM_ITER_ELEM (f, &fiter, e, BM_FACES_OF_EDGE) {
      /* could support ngons, other areas would need updating too, see T48926. */
      if ((f->len <= 4) && !BMO_face_flag_test(bm, f, FACE_OUT)) {
        BMIter liter;
        BMLoop *l;
        bool ok = false;
 
        /* check at least 2 edges in the face are rings */
        BM_ITER_ELEM (l, &liter, f, BM_LOOPS_OF_FACE) {
          if (BMO_edge_flag_test(bm, l->e, EDGE_RING) && e != l->e) {
            ok = true;
            break;
          }
        }
 
        if (ok) {
          BMO_face_flag_enable(bm, f, FACE_OUT);
 
          BM_ITER_ELEM (l, &liter, f, BM_LOOPS_OF_FACE) {
            if (!BMO_edge_flag_test(bm, l->e, EDGE_RING)) {
              BMO_edge_flag_enable(bm, l->e, EDGE_RIM);
            }
          }
        }
      }
    }
  }
 
  /* -------------------------------------------------------------------- */
  /* Cache falloff for each step (symmetrical) */
 
  if (falloff_cache) {
    int i;
    for (i = 0; i < resolu; i++) {
      float shape_size = 1.0f;
      float fac = (float)i / (float)(resolu - 1);
      fac = fabsf(1.0f - 2.0f * fabsf(0.5f - fac));
      fac = bmesh_subd_falloff_calc(profile_shape, fac);
      shape_size += fac * profile_shape_factor;
 
      falloff_cache[i] = shape_size;
    }
  }
 
  /* -------------------------------------------------------------------- */
  /* Execute subdivision on all ring pairs */
 
  count = BM_mesh_edgeloops_find(bm, &eloops_rim, bm_edge_rim_test_cb, (void *)bm);
 
  if (count < 2) {
    BMO_error_raise(bm, op, BMERR_INVALID_SELECTION, "No edge rings found");
    goto cleanup;
  }
  else if (count == 2) {
    /* this case could be removed,
     * but simple to avoid 'bm_edgering_pair_calc' in this case since there's only one. */
    struct BMEdgeLoopStore *el_store_a = eloops_rim.first;
    struct BMEdgeLoopStore *el_store_b = eloops_rim.last;
    LoopPairStore *lpair;
 
    if (bm_edgeloop_check_overlap_all(bm, el_store_a, el_store_b)) {
      lpair = bm_edgering_pair_store_create(bm, el_store_a, el_store_b, interp_mode);
    }
    else {
      lpair = NULL;
    }
 
    if (lpair) {
      bm_edgering_pair_ringsubd(
          bm, lpair, el_store_a, el_store_b, interp_mode, cuts, smooth, falloff_cache);
      bm_edgering_pair_store_free(lpair, interp_mode);
      changed = true;
    }
    else {
      BMO_error_raise(bm, op, BMERR_INVALID_SELECTION, "Edge-ring pair isn't connected");
      goto cleanup;
    }
  }
  else {
    GSetIterator gs_iter;
    int i;
 
    GSet *eloop_pairs_gs = bm_edgering_pair_calc(bm, &eloops_rim);
    LoopPairStore **lpair_arr;
 
    if (eloop_pairs_gs == NULL) {
      BMO_error_raise(bm, op, BMERR_INVALID_SELECTION, "Edge-rings are not connected");
      goto cleanup;
    }
 
    lpair_arr = BLI_array_alloca(lpair_arr, BLI_gset_len(eloop_pairs_gs));
 
    /* first cache pairs */
    GSET_ITER_INDEX (gs_iter, eloop_pairs_gs, i) {
      GHashPair *eloop_pair = BLI_gsetIterator_getKey(&gs_iter);
      struct BMEdgeLoopStore *el_store_a = (void *)eloop_pair->first;
      struct BMEdgeLoopStore *el_store_b = (void *)eloop_pair->second;
      LoopPairStore *lpair;
 
      if (bm_edgeloop_check_overlap_all(bm, el_store_a, el_store_b)) {
        lpair = bm_edgering_pair_store_create(bm, el_store_a, el_store_b, interp_mode);
      }
      else {
        lpair = NULL;
      }
      lpair_arr[i] = lpair;
 
      BLI_assert(bm_verts_tag_count(bm) == 0);
    }
 
    GSET_ITER_INDEX (gs_iter, eloop_pairs_gs, i) {
      GHashPair *eloop_pair = BLI_gsetIterator_getKey(&gs_iter);
      struct BMEdgeLoopStore *el_store_a = (void *)eloop_pair->first;
      struct BMEdgeLoopStore *el_store_b = (void *)eloop_pair->second;
      LoopPairStore *lpair = lpair_arr[i];
 
      if (lpair) {
        bm_edgering_pair_ringsubd(
            bm, lpair, el_store_a, el_store_b, interp_mode, cuts, smooth, falloff_cache);
        bm_edgering_pair_store_free(lpair, interp_mode);
        changed = true;
      }
 
      BLI_assert(bm_verts_tag_count(bm) == 0);
    }
    BLI_gset_free(eloop_pairs_gs, MEM_freeN);
  }
 
cleanup:
  BM_mesh_edgeloops_free(&eloops_rim);
 
  /* flag output */
  if (changed) {
    BMO_slot_buffer_from_enabled_flag(bm, op, op->slots_out, "faces.out", BM_FACE, FACE_OUT);
  }
}