bmo_join_triangles.cc

Go to the documentation of this file.

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

 
#include <algorithm>
 
#include "MEM_guardedalloc.h"
 
#include "BLI_heap.h"
#include "BLI_math_base.h"
#include "BLI_math_geom.h"
#include "BLI_math_rotation.h"
#include "BLI_math_vector.h"
 
#include "BKE_customdata.hh"
 
#include "bmesh.hh"
 
#include "intern/bmesh_operators_private.hh" /* own include */

#define ASSERT_VALID_ERROR_METRIC(val) BLI_assert(isfinite(val) && (val) >= 0 && (val) <= 2 * M_PI)
 
#if 0
#  define USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
#endif
 
#define FACE_OUT (1 << 0)
#define FACE_INPUT (1 << 2)

constexpr float maximum_improvement = 0.99f;
 
/* -------------------------------------------------------------------- */
 
struct JoinEdgesState {
  Heap *edge_queue;

  HeapNode **edge_queue_nodes;

  bool use_topo_influence;

  float topo_influnce;

  bool select_tris_only;
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  BMesh *debug_bm;

  int debug_merge_limit;

  int debug_merge_count;

  int debug_neighbor;

  int debug_neighbor_global_count;

  bool debug_this_step;
#endif
};

/* -------------------------------------------------------------------- */

static float quad_calc_error(const float v1[3],
                             const float v2[3],
                             const float v3[3],
                             const float v4[3])
{
  float error = 0.0f;
 
  /* Normal difference: a perfectly flat planar face adds a difference of zero. */
  {
    float n1[3], n2[3];
    float angle_a, angle_b;
    float diff;
 
    normal_tri_v3(n1, v1, v2, v3);
    normal_tri_v3(n2, v1, v3, v4);
    angle_a = compare_v3v3(n1, n2, FLT_EPSILON) ? 0.0f : angle_normalized_v3v3(n1, n2);
 
    normal_tri_v3(n1, v2, v3, v4);
    normal_tri_v3(n2, v4, v1, v2);
    angle_b = compare_v3v3(n1, n2, FLT_EPSILON) ? 0.0f : angle_normalized_v3v3(n1, n2);
 
    diff = (angle_a + angle_b) / float(M_PI * 2);
 
    ASSERT_VALID_ERROR_METRIC(diff);
 
    error += diff;
  }
 
  /* Co-linearity: a face with four right angle corners adds a difference of zero. */
  {
    float edge_vecs[4][3];
    float diff;
 
    sub_v3_v3v3(edge_vecs[0], v1, v2);
    sub_v3_v3v3(edge_vecs[1], v2, v3);
    sub_v3_v3v3(edge_vecs[2], v3, v4);
    sub_v3_v3v3(edge_vecs[3], v4, v1);
 
    normalize_v3(edge_vecs[0]);
    normalize_v3(edge_vecs[1]);
    normalize_v3(edge_vecs[2]);
    normalize_v3(edge_vecs[3]);
 
    /* A completely skinny face is 'pi' after halving. */
    diff = (fabsf(angle_normalized_v3v3(edge_vecs[0], edge_vecs[1]) - float(M_PI_2)) +
            fabsf(angle_normalized_v3v3(edge_vecs[1], edge_vecs[2]) - float(M_PI_2)) +
            fabsf(angle_normalized_v3v3(edge_vecs[2], edge_vecs[3]) - float(M_PI_2)) +
            fabsf(angle_normalized_v3v3(edge_vecs[3], edge_vecs[0]) - float(M_PI_2))) /
           float(M_PI * 2);
 
    ASSERT_VALID_ERROR_METRIC(diff);
 
    error += diff;
  }
 
  /* Concavity: a face with no concavity adds an error of 0. */
  {
    float area_min, area_max, area_a, area_b;
    float diff;
 
    area_a = area_tri_v3(v1, v2, v3) + area_tri_v3(v1, v3, v4);
    area_b = area_tri_v3(v2, v3, v4) + area_tri_v3(v4, v1, v2);
 
    area_min = min_ff(area_a, area_b);
    area_max = max_ff(area_a, area_b);
 
    /* Note use of ternary operator to guard against divide by zero. */
    diff = area_max ? (1.0f - (area_min / area_max)) : 1.0f;
 
    ASSERT_VALID_ERROR_METRIC(diff);
 
    error += diff;
  }
 
  ASSERT_VALID_ERROR_METRIC(error);
 
  return error;
}

static void bm_edge_to_quad_verts(const BMEdge *e, const BMVert *r_v_quad[4])
{
  BLI_assert(BM_edge_is_manifold(e));
  BLI_assert(e->l->f->len == 3 && e->l->radial_next->f->len == 3);
  r_v_quad[0] = e->l->v;
  r_v_quad[1] = e->l->radial_next->prev->v;
  r_v_quad[2] = e->l->next->v;
  r_v_quad[3] = e->l->prev->v;
}
 
/* -------------------------------------------------------------------- */

 
namespace {
 
struct DelimitData_CD {
  int cd_type;
  int cd_size;
  int cd_offset;
  int cd_offset_end;
};
 
struct DelimitData {
  uint do_seam : 1;
  uint do_sharp : 1;
  uint do_mat : 1;
  uint do_angle_face : 1;
  uint do_angle_shape : 1;
 
  float angle_face;
  float angle_face__cos;
 
  float angle_shape;
 
  DelimitData_CD cdata[4];
  int cdata_len;
};
 
}  // namespace

static bool bm_edge_is_contiguous_loop_cd_all(const BMEdge *e, const DelimitData_CD *delimit_data)
{
  int cd_offset;
  for (cd_offset = delimit_data->cd_offset; cd_offset < delimit_data->cd_offset_end;
       cd_offset += delimit_data->cd_size)
  {
    if (BM_edge_is_contiguous_loop_cd(e, delimit_data->cd_type, cd_offset) == false) {
      return false;
    }
  }
 
  return true;
}

static bool bm_edge_delimit_cdata(CustomData *ldata,
                                  eCustomDataType type,
                                  DelimitData_CD *r_delim_cd)
{
  const int layer_len = CustomData_number_of_layers(ldata, type);
  r_delim_cd->cd_type = type;
  r_delim_cd->cd_size = CustomData_sizeof(eCustomDataType(r_delim_cd->cd_type));
  r_delim_cd->cd_offset = CustomData_get_n_offset(ldata, type, 0);
  r_delim_cd->cd_offset_end = r_delim_cd->cd_offset + (r_delim_cd->cd_size * layer_len);
  return (r_delim_cd->cd_offset != -1);
}

static DelimitData bm_edge_delmimit_data_from_op(BMesh *bm, BMOperator *op)
{
  DelimitData delimit_data = {0};
  delimit_data.do_seam = BMO_slot_bool_get(op->slots_in, "cmp_seam");
  delimit_data.do_sharp = BMO_slot_bool_get(op->slots_in, "cmp_sharp");
  delimit_data.do_mat = BMO_slot_bool_get(op->slots_in, "cmp_materials");
 
  /* Determine if angle face processing occurs and its parameters. */
  float angle_face = BMO_slot_float_get(op->slots_in, "angle_face_threshold");
  if (angle_face < DEG2RADF(180.0f)) {
    delimit_data.angle_face = angle_face;
    delimit_data.angle_face__cos = cosf(angle_face);
    delimit_data.do_angle_face = true;
  }
  else {
    delimit_data.do_angle_face = false;
  }
 
  /* Determine if angle shape processing occurs and its parameters. */
  float angle_shape = BMO_slot_float_get(op->slots_in, "angle_shape_threshold");
  if (angle_shape < DEG2RADF(180.0f)) {
    delimit_data.angle_shape = angle_shape;
    delimit_data.do_angle_shape = true;
  }
  else {
    delimit_data.do_angle_shape = false;
  }
 
  if (BMO_slot_bool_get(op->slots_in, "cmp_uvs") &&
      bm_edge_delimit_cdata(
          &bm->ldata, CD_PROP_FLOAT2, &delimit_data.cdata[delimit_data.cdata_len]))
  {
    delimit_data.cdata_len += 1;
  }
 
  delimit_data.cdata[delimit_data.cdata_len].cd_offset = -1;
  if (BMO_slot_bool_get(op->slots_in, "cmp_vcols") &&
      bm_edge_delimit_cdata(
          &bm->ldata, CD_PROP_BYTE_COLOR, &delimit_data.cdata[delimit_data.cdata_len]))
  {
    delimit_data.cdata_len += 1;
  }
  return delimit_data;
}

static bool bm_edge_is_delimit(const BMEdge *e, const DelimitData *delimit_data)
{
  BMFace *f_a = e->l->f, *f_b = e->l->radial_next->f;
#if 0
  const bool is_contig = BM_edge_is_contiguous(e);
  float angle;
#endif
 
  if (delimit_data->do_seam && BM_elem_flag_test(e, BM_ELEM_SEAM)) {
    return true;
  }
 
  if (delimit_data->do_sharp && (BM_elem_flag_test(e, BM_ELEM_SMOOTH) == 0)) {
    return true;
  }
 
  if (delimit_data->do_mat && (f_a->mat_nr != f_b->mat_nr)) {
    return true;
  }
 
  if (delimit_data->do_angle_face) {
    if (dot_v3v3(f_a->no, f_b->no) < delimit_data->angle_face__cos) {
      return true;
    }
  }
 
  if (delimit_data->do_angle_shape) {
    const BMVert *verts[4];
    bm_edge_to_quad_verts(e, verts);
 
    /* if we're checking the shape at all, a flipped face is out of the question */
    if (is_quad_flip_v3(verts[0]->co, verts[1]->co, verts[2]->co, verts[3]->co)) {
      return true;
    }
 
    float edge_vecs[4][3];
 
    sub_v3_v3v3(edge_vecs[0], verts[0]->co, verts[1]->co);
    sub_v3_v3v3(edge_vecs[1], verts[1]->co, verts[2]->co);
    sub_v3_v3v3(edge_vecs[2], verts[2]->co, verts[3]->co);
    sub_v3_v3v3(edge_vecs[3], verts[3]->co, verts[0]->co);
 
    normalize_v3(edge_vecs[0]);
    normalize_v3(edge_vecs[1]);
    normalize_v3(edge_vecs[2]);
    normalize_v3(edge_vecs[3]);
 
    if ((fabsf(angle_normalized_v3v3(edge_vecs[0], edge_vecs[1]) - float(M_PI_2)) >
         delimit_data->angle_shape) ||
        (fabsf(angle_normalized_v3v3(edge_vecs[1], edge_vecs[2]) - float(M_PI_2)) >
         delimit_data->angle_shape) ||
        (fabsf(angle_normalized_v3v3(edge_vecs[2], edge_vecs[3]) - float(M_PI_2)) >
         delimit_data->angle_shape) ||
        (fabsf(angle_normalized_v3v3(edge_vecs[3], edge_vecs[0]) - float(M_PI_2)) >
         delimit_data->angle_shape))
    {
      return true;
    }
  }
 
  if (delimit_data->cdata_len) {
    int i;
    for (i = 0; i < delimit_data->cdata_len; i++) {
      if (!bm_edge_is_contiguous_loop_cd_all(e, &delimit_data->cdata[i])) {
        return true;
      }
    }
  }
 
  return false;
}

 
struct JoinEdgesNeighborItem {
  BMEdge *e;
  BMLoop *l;
};
 
struct JoinEdgesNeighborInfo {
  JoinEdgesNeighborItem items[8];
  int items_num;
};

static void add_without_duplicates(JoinEdgesNeighborInfo &neighbor_info, BMEdge *e, BMLoop *l)
{
  BLI_assert(neighbor_info.items_num < ARRAY_SIZE(neighbor_info.items));
 
  /* Don't add null pointers. Another safeguard which cannot happen. */
  BLI_assert(e != nullptr);
 
  /* Don't add duplicates. */
  for (uint index = 0; index < neighbor_info.items_num; index++) {
    if (neighbor_info.items[neighbor_info.items_num].e == e) {
      return;
    }
  }
 
  /* Add the edge and increase the count by 1. */
  JoinEdgesNeighborItem *item = &neighbor_info.items[neighbor_info.items_num++];
  item->e = e;
  item->l = l;
}

static void add_neighbors(JoinEdgesNeighborInfo &neighbor_info, BMLoop *l_in_quad)
{
  /* If the edge is not manifold, there is no neighboring face to process. */
  if (!BM_edge_is_manifold(l_in_quad->e)) {
    /* No new edges added. */
    return;
  }
 
  BMLoop *l_in_neighbor = l_in_quad->radial_next;
 
  /* If the neighboring face is not a triangle, don't process it. */
  if (l_in_neighbor->f->len != 3) {
    /* No new edges added. */
    return;
  }
 
#ifndef NDEBUG
  const int items_num_prev = neighbor_info.items_num;
#endif
 
  /* Get the other two loops of the neighboring triangle. */
  BMLoop *l_other_arr[2] = {l_in_neighbor->prev, l_in_neighbor->next};
  for (int i = 0; i < ARRAY_SIZE(l_other_arr); i++) {
    BMLoop *l_other = l_other_arr[i];
 
    /* If `l_other` is manifold, and the adjacent face is also a triangle,
     * mark it for potential improvement. */
    if (BM_edge_is_manifold(l_other->e) && l_other->radial_next->f->len == 3) {
      add_without_duplicates(neighbor_info, l_other->e, l_in_neighbor);
    }
  }
 
  /* Added either 0, 1, or 2 edges. */
#ifndef NDEBUG
  BLI_assert(neighbor_info.items_num - items_num_prev < 3);
#endif
}

static void rotate_to_plane(const JoinEdgesState &s,
                            const BMVert *quad_verts[4],
                            const BMLoop *l_shared,
                            const float plane_normal[3],
                            float r_quad_coordinates[4][3])
{
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  if (s.debug_this_step) {
    printf("angle");
  }
#else
  UNUSED_VARS(s);
#endif
 
  float rotation_axis[3];
  sub_v3_v3v3(rotation_axis, l_shared->v->co, l_shared->next->v->co);
  normalize_v3(rotation_axis);
 
  float quad_normal[3] = {0};
  normal_quad_v3(
      quad_normal, quad_verts[0]->co, quad_verts[1]->co, quad_verts[2]->co, quad_verts[3]->co);
 
  float angle = angle_signed_on_axis_v3v3_v3(plane_normal, quad_normal, rotation_axis);
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  if (s.debug_this_step) {
    printf(" %f, ", RAD2DEGF(angle));
  }
#endif
 
  for (int i = 0; i < 4; i++) {
    if (ELEM(quad_verts[i], l_shared->v, l_shared->next->v)) {
      /* Two coordinates of the quad match the vector that defines the axis of rotation, so they
       * don't change. */
      copy_v3_v3(r_quad_coordinates[i], quad_verts[i]->co);
    }
    else {
      /* The other two coordinates get rotated around the axis, and so they change. */
      float local_coordinate[3];
      sub_v3_v3v3(local_coordinate, quad_verts[i]->co, l_shared->v->co);
      rotate_normalized_v3_v3v3fl(r_quad_coordinates[i], local_coordinate, rotation_axis, angle);
      add_v3_v3(r_quad_coordinates[i], l_shared->v->co);
    }
  }
}

static float compute_alignment(const JoinEdgesState &s,
                               const float quad_a_vecs[4][3],
                               const BMVert *quad_b_verts[4],
                               const BMLoop *l_shared,
                               const float plane_normal[3])
{
  /* Many meshes have lots of curvature or sharp edges. Pairs of quads shouldn't be penalized
   * *worse* because they represent a curved surface or define an edge. So we rotate quad_b around
   * its common edge with quad_a until both are, as much as possible, in the same plane.
   * This ensures the best possible chance to align. */
  float quad_b_coordinates[4][3];
  rotate_to_plane(s, quad_b_verts, l_shared, plane_normal, quad_b_coordinates);
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  if (s.debug_this_step) {
    /* For visualization purposes *only*, rotate the face being considered.
     * The const_cast here is purposeful. We want to specify `const BMVert` deliberately
     * to show that we're not *supposed* to be moving verts around.
     * But only for debug visualization, we do. This alters the mesh to visualize the
     * effect of rotating the face into the plane for alignment testing. */
    copy_v3_v3(const_cast<BMVert *>(quad_b_verts[0])->co, quad_b_coordinates[0]);
    copy_v3_v3(const_cast<BMVert *>(quad_b_verts[1])->co, quad_b_coordinates[1]);
    copy_v3_v3(const_cast<BMVert *>(quad_b_verts[2])->co, quad_b_coordinates[2]);
    copy_v3_v3(const_cast<BMVert *>(quad_b_verts[3])->co, quad_b_coordinates[3]);
  }
#endif
 
  /* compute the four unit vectors of the quad b edges. */
  float quad_b_vecs[4][3];
  sub_v3_v3v3(quad_b_vecs[0], quad_b_coordinates[0], quad_b_coordinates[1]);
  sub_v3_v3v3(quad_b_vecs[1], quad_b_coordinates[1], quad_b_coordinates[2]);
  sub_v3_v3v3(quad_b_vecs[2], quad_b_coordinates[2], quad_b_coordinates[3]);
  sub_v3_v3v3(quad_b_vecs[3], quad_b_coordinates[3], quad_b_coordinates[0]);
  normalize_v3(quad_b_vecs[0]);
  normalize_v3(quad_b_vecs[1]);
  normalize_v3(quad_b_vecs[2]);
  normalize_v3(quad_b_vecs[3]);
 
  /* Given that we're not certain of how the first loop of the quad and the first loop
   * of the proposed merge quad relate to each other, there are four possible combinations
   * to check, to test that the neighbor face and the merged face have good alignment.
   *
   * In theory, a very nuanced analysis involving l_shared, loop pointers, vertex pointers,
   * etc, would allow determining which sets of vectors are the right matches sets to compare.
   *
   * Do not meddle in the affairs of algorithms, for they are subtle and quick to anger.
   *
   * Instead, this code does the math twice, then it just flips each component by 180 degrees to
   * pick up the other two cases. Four extra angle tests aren't that much worse than optimal.
   * Brute forcing the math and ending up with clear and understandable code is better. */
 
  float error[4] = {0.0f};
  for (int i = 0; i < ARRAY_SIZE(error); i++) {
    const float angle_a = fabsf(angle_normalized_v3v3(quad_a_vecs[i], quad_b_vecs[i]));
    const float angle_b = fabsf(angle_normalized_v3v3(quad_a_vecs[i], quad_b_vecs[(i + 1) % 4]));
 
    /* Compute the case if the quads are aligned. */
    error[0] += angle_a;
    /* Compute the case if the quads are 90 degrees rotated. */
    error[1] += angle_b;
 
    /* Compute the case if the quads are 180 degrees rotated.
     * This is `error[0]` except each error component is individually rotated 180 degrees. */
    error[2] += M_PI - angle_a;
 
    /* Compute the case if the quads are 270 degrees rotated.
     * This is `error[1]` except each error component is individually rotated 180 degrees. */
    error[3] += M_PI - angle_b;
  }
 
  /* Pick the best option and average the four components. */
  const float best_error = std::min({error[0], error[1], error[2], error[3]}) / 4.0f;
 
  ASSERT_VALID_ERROR_METRIC(best_error);
 
  /* Based on the best error, we scale how aligned we are to the range 0...1
   * `M_PI / 4` is used here because the worst case is a quad with all four edges
   * at 45 degree angles. */
  float alignment = 1.0f - (best_error / (M_PI / 4.0f));
 
  /* if alignment is *truly* awful, then do nothing. Don't make a join worse. */
  alignment = std::max(alignment, 0.0f);
 
  ASSERT_VALID_ERROR_METRIC(alignment);
 
  return alignment;
}

static void reprioritize_join(JoinEdgesState &s,
                              BMEdge *e_merge,
                              BMLoop *l_shared,
                              float neighbor_quad_vecs[4][3],
                              const float neighbor_quad_error,
                              const float neighbor_quad_normal[3])
{
  ASSERT_VALID_ERROR_METRIC(neighbor_quad_error);
 
  /* If the edge wasn't found, (delimit, non-manifold, etc) then return.
   * Nothing to do here. */
  BLI_assert(BM_elem_index_get(e_merge) >= 0);
  HeapNode *node = s.edge_queue_nodes[BM_elem_index_get(e_merge)];
  if (node == nullptr) {
    return;
  }
 
  float join_error_curr = BLI_heap_node_value(node);
 
  ASSERT_VALID_ERROR_METRIC(join_error_curr);
 
  /* Never make a join *worse* because of topology around it.
   * Because we are sorted during the join phase of the algorithm, this should *only* happen when
   * processing any pre-existing quads in the input mesh during setup. They might have high error.
   * If they do, ignore them. */
  if (neighbor_quad_error > join_error_curr) {
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
    /* This should only happen during setup.
     * Indicates an error, if it happens once we've started merging. */
    BLI_assert(s.debug_merge_count == 0);
#endif
 
    return;
  }
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  if (s.debug_this_step) {
    printf("Edge improved from ");
  }
#endif
 
  /* Get the four corners of the quad that would result if we merged. */
  const BMVert *quad_verts_merge[4];
  bm_edge_to_quad_verts(e_merge, quad_verts_merge);
 
  /* Now compute the alignment.
   * Regular grids of rectangles or trapezoids have high alignment
   * Mismatched combinations of rectangles diamonds and trapezoids have low alignment. */
  const float alignment = compute_alignment(
      s, neighbor_quad_vecs, quad_verts_merge, l_shared, neighbor_quad_normal);
 
  /* Compute how much the neighbor is better than the candidate.
   * Since the neighbor quad error is smaller, Improvement is always represented as negative. */
  const float improvement = (neighbor_quad_error - join_error_curr);
 
  ASSERT_VALID_ERROR_METRIC(-improvement);
 
  /* Compute the scale factor for how much of that possible improvement we should apply to this
   * edge. This combines... topology_influence, which is an operator setting... and alignment,
   * which is computed. Faces which are diagonal have an alignment of 0% - perfect rectangular
   * grids have an alignment of 100% Neither topology_influence nor alignment can be negative;
   * therefore the multiplier *never* makes error worse. once combined, 0 means no improvement, 1
   * means improve all the way to exactly match the quality of the contributing neighbor.
   * topology_influece is allowed to exceed 1.0, which lets it cancel out some of the alignment
   * penalty. */
  float multiplier = s.topo_influnce * alignment;
 
  /* However, the combined multiplier shouldn't ever be allowed to exceed 1.0 because permitting
   * that would cause exponential growth when alignment is very good, and when that happens, the
   * algorithm becomes crazy.
   *
   * Further, if we allow a multiplier of exactly 1.0, then all eight edges around the neighbor
   * quad would end up with a quality that is *exactly* equal to the neighbor - and each other;
   * losing valuable information about their relative sorting.
   * In order to preserve that, the multiplier is capped at 99%.
   * The last 1% that is left uncorrected is enough to preserve relative ordering.
   *
   * This especially helps in quads that touch 3-poles and 5-poles. Since those quads naturally
   * have diamond shapes, their initial error values tend to be higher and they sort to the end of
   * the priority queue. Limiting improvement at 99% ensures those quads tend to retain their bad
   * sort, meaning they end up surrounded by quads that define a good grid,
   * then they merge last, which tends to produce better results. */
  multiplier = std::min(multiplier, maximum_improvement);
 
  ASSERT_VALID_ERROR_METRIC(multiplier);
 
  /* improvement is always represented as a negative number (that will reduce error)
   * Based on that convention, `+` is correct here. */
  float join_error_next = join_error_curr + (improvement * multiplier);
 
  ASSERT_VALID_ERROR_METRIC(join_error_next);
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  if (s.debug_this_step) {
    printf("%f to %f with alignment of %f.\n", join_error_curr, join_error_next, alignment);
    BMO_face_flag_enable(s.debug_bm, e_merge->l->f, FACE_OUT);
    BMO_face_flag_enable(s.debug_bm, e_merge->l->radial_next->f, FACE_OUT);
  }
#endif
 
  /* Now, update the node value in the heap, which may cause the node
   * to be moved toward the head of the priority queue. */
  BLI_heap_node_value_update(s.edge_queue, node, join_error_next);
}

static void reprioritize_face_neighbors(JoinEdgesState &s, BMFace *f, float f_error)
{
  BLI_assert(f->len == 4);
 
  /* Identify any mergeable edges of any neighbor triangles that face us.
   * - Some of our four edges... might not be manifold.
   * - Some of our neighbor faces... might not be triangles.
   * - Some of our neighbor triangles... might have other non-manifold (unmergeable) edges.
   * - Some of our neighbor triangles' manifold edges... might have non-triangle neighbors.
   * Therefore, there can be have up to eight mergeable edges, although there are often fewer. */
  JoinEdgesNeighborInfo neighbor_info = {};
 
  /* Get the four loops around the face. */
  BMLoop *l_quad[4];
  BM_face_as_array_loop_quad(f, l_quad);
 
  /* Add the mergeable neighbors for each of those loops. */
  for (int i = 0; i < ARRAY_SIZE(l_quad); i++) {
    add_neighbors(neighbor_info, l_quad[i]);
  }
 
  /* Return if there is nothing to do. */
  if (neighbor_info.items_num == 0) {
    return;
  }
 
  /* Compute the four unit vectors around this quad. */
  float quad_vecs[4][3];
  for (int i_next = 0, i = ARRAY_SIZE(l_quad) - 1; i_next < ARRAY_SIZE(l_quad); i = i_next++) {
    sub_v3_v3v3(quad_vecs[i], l_quad[i]->v->co, l_quad[i_next]->v->co);
    normalize_v3(quad_vecs[i]);
  }
 
  /* Re-prioritize each neighbor. */
  for (int i = 0; i < neighbor_info.items_num; i++) {
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
    s.debug_this_step = (s.debug_merge_limit > 0 && s.debug_merge_count == s.debug_merge_limit &&
                         i + 1 == s.debug_neighbor) ||
                        (++s.debug_neighbor_global_count == s.debug_neighbor &&
                         s.debug_merge_limit == 0);
#endif
    const JoinEdgesNeighborItem *item = &neighbor_info.items[i];
    reprioritize_join(s, item->e, item->l, quad_vecs, f_error, f->no);
  }
}

static BMFace *bm_faces_join_pair_by_edge(BMesh *bm,
                                          BMEdge *e
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
                                          ,
                                          JoinEdgesState &s
#endif
)
{
  /* Non-manifold edges can't be merged. */
  BLI_assert(BM_edge_is_manifold(e));
 
  /* Identify the loops on either side of the edge which may be joined. */
  BMLoop *l_a = e->l;
  BMLoop *l_b = e->l->radial_next;
 
  /* If previous face merges have created quads, which now make this edge unmergeable,
   * then skip it and move on. This happens frequently and that's ok.
   * It's much easier and more efficient to just skip these edges when we encounter them,
   * than it is to try to search the heap for them and remove them preemptively. */
  if ((l_a->f->len != 3) || (l_b->f->len != 3)) {
    return nullptr;
  }
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  /* Stop doing merges in the middle of processing if we reached a user limit.
   * This is allowed so a developer can check steps in the process of the algorithm. */
  if (++s.debug_merge_count > s.debug_merge_limit && s.debug_merge_limit != -1) {
    return nullptr;
  }
#endif
 
  BMFace *f_double;
 
  /* Join the edge and identify the face. */
  BMFace *f = BM_faces_join_pair(bm, l_a, l_b, true, &f_double);
  /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
  BLI_assert_msg(f_double == nullptr,
                 "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
 
  return f;
}

void bmo_join_triangles_exec(BMesh *bm, BMOperator *op)
{
  BMIter iter;
  BMOIter siter;
  BMFace *f;
  BMEdge *e;
 
  DelimitData delimit_data = bm_edge_delmimit_data_from_op(bm, op);
 
  /* Initial setup of state. */
  JoinEdgesState s = {nullptr};
  s.topo_influnce = BMO_slot_float_get(op->slots_in, "topology_influence");
  s.use_topo_influence = (s.topo_influnce != 0.0f);
  s.edge_queue = BLI_heap_new();
  s.select_tris_only = BMO_slot_bool_get(op->slots_in, "deselect_joined");
  if (s.use_topo_influence) {
    s.edge_queue_nodes = MEM_malloc_arrayN<HeapNode *>(bm->totedge, __func__);
  }
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  s.debug_bm = bm;
  s.debug_merge_limit = BMO_slot_int_get(op->slots_in, "merge_limit") - 1;
  s.debug_neighbor = BMO_slot_int_get(op->slots_in, "neighbor_debug");
  s.debug_merge_count = 0;
  s.debug_neighbor_global_count = 0;
#endif
 
  /* Go through every face in the input slot. Mark triangles for processing. */
  BMO_ITER (f, &siter, op->slots_in, "faces", BM_FACE) {
    if (f->len == 3) {
      BMO_face_flag_enable(bm, f, FACE_INPUT);
 
      /* And setup the initial selection. */
      if (s.select_tris_only) {
        BMO_face_flag_enable(bm, f, FACE_OUT);
      }
    }
  }
 
  /* Go through every edge in the mesh, mark edges that can be merged. */
  int i = 0;
  BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
    BM_elem_index_set(e, i); /* set_inline */
 
    /* If the edge is manifold, has a tagged input triangle on both sides,
     * and is *not* delimited, then it's a candidate to merge. */
    BMFace *f_a, *f_b;
    if (BM_edge_face_pair(e, &f_a, &f_b) && BMO_face_flag_test(bm, f_a, FACE_INPUT) &&
        BMO_face_flag_test(bm, f_b, FACE_INPUT) && !bm_edge_is_delimit(e, &delimit_data))
    {
      /* Compute the error that would result from a merge. */
      const BMVert *e_verts[4];
      bm_edge_to_quad_verts(e, e_verts);
      const float merge_error = quad_calc_error(
          e_verts[0]->co, e_verts[1]->co, e_verts[2]->co, e_verts[3]->co);
 
      /* Record the candidate merge in both the heap, and the heap index. */
      HeapNode *node = BLI_heap_insert(s.edge_queue, merge_error, e);
      if (s.use_topo_influence) {
        s.edge_queue_nodes[i] = node;
      }
    }
    else {
      if (s.use_topo_influence) {
        s.edge_queue_nodes[i] = nullptr;
      }
    }
  }
 
  /* Go through all the faces of the input slot, this time to find quads.
   * Improve the candidates around any preexisting quads in the mesh.
   *
   * NOTE: This unfortunately misses any quads which are not selected, but
   * which neighbor the selection. The only alternate would be to iterate the
   * whole mesh, which might be expensive for very large meshes with small selections. */
  if (s.use_topo_influence && (BLI_heap_is_empty(s.edge_queue) == false)) {
    BMO_ITER (f, &siter, op->slots_in, "faces", BM_FACE) {
      if (f->len == 4) {
        BMVert *f_verts[4];
        BM_face_as_array_vert_quad(f, f_verts);
 
        /* Flat quads with right angle corners and no concavity have lower error. */
        float f_error = quad_calc_error(
            f_verts[0]->co, f_verts[1]->co, f_verts[2]->co, f_verts[3]->co);
 
        /* Apply the compensated error.
         * Since we're early in the process we over-prioritize any already existing quads to
         * allow them to have an especially strong influence on the resulting mesh.
         * At a topology influence of 200%, they're considered to be *almost perfect* quads
         * regardless of their actual error. Either way, the multiplier is never completely
         * allowed to reach zero. Instead, 1% of the original error is preserved...
         * which is enough to maintain the relative priority sorting between existing quads. */
        f_error *= (2.0f - (s.topo_influnce * maximum_improvement));
 
        reprioritize_face_neighbors(s, f, f_error);
      }
    }
  }
 
  /* Process all possible merges. */
  while (!BLI_heap_is_empty(s.edge_queue)) {
 
    /* Get the best merge from the priority queue.
     * Remove it from the priority queue. */
    const float f_error = BLI_heap_top_value(s.edge_queue);
    BMEdge *e = reinterpret_cast<BMEdge *>(BLI_heap_pop_min(s.edge_queue));
 
    /* Attempt the merge. */
    BMFace *f_new = bm_faces_join_pair_by_edge(bm,
                                               e
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
                                               ,
                                               s
#endif
    );
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
    /* If stopping partway through, clear the selection entirely, and instead
     * highlight the faces being considered in the step the user is checking. */
    if (s.debug_merge_limit != -1 && s.debug_merge_count == s.debug_merge_limit) {
      BMEdge *f;
      BMIter iter;
      BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
        BMO_face_flag_disable(bm, f, FACE_OUT);
      }
    }
#endif
 
    if (f_new) {
 
      /* Tag the face so the selection can be extended to include the new face. */
      if (s.select_tris_only == false) {
        BMO_face_flag_enable(bm, f_new, FACE_OUT);
      }
 
      /* Improve the neighbors on success. */
      if (s.use_topo_influence) {
        reprioritize_face_neighbors(s, f_new, f_error);
      }
    }
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
    /* If we're supposed to stop partway though, do that.
     * This allows a developer to inspect the mesh at intermediate stages of processing. */
    if (s.debug_merge_limit != -1 && s.debug_merge_count >= s.debug_merge_limit) {
      break;
    }
#endif
  }
 
#ifdef USE_JOIN_TRIANGLE_INTERACTIVE_TESTING
  /* Expect a full processing to have occurred, *only* if we didn't stop partway through. */
  if (!(s.debug_merge_limit != -1 && s.debug_merge_count >= s.debug_merge_limit))
#endif
  {
    /* Expect a full processing to have occurred. */
    BLI_assert(BLI_heap_is_empty(s.edge_queue));
  }
 
  /* Clean up. */
  BLI_heap_free(s.edge_queue, nullptr);
  if (s.use_topo_influence) {
    MEM_freeN(s.edge_queue_nodes);
  }
 
  /* Return the selection results. */
  BMO_slot_buffer_from_enabled_flag(bm, op, op->slots_out, "faces.out", BM_FACE, FACE_OUT);
}