text: /*
   * 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.
   *
   * The Original Code is Copyright (C) 2019 Blender Foundation.
   * All rights reserved.
   */
  
  /** \file
   * \ingroup bmesh
   */
  
  #include "MEM_guardedalloc.h"
  
  #include "BLI_math.h"
  #include "BLI_sort.h"
  #include "BLI_stack.h"
  
  #include "BKE_bvhutils.h"
  
  #include "bmesh.h"
  
  #include "bmesh_intersect_edges.h" /* own include */
  
  #define KDOP_AXIS_LEN 14
  
  /* -------------------------------------------------------------------- */
  /** \name Weld Linked Wire Edges into Linked Faces
   *
   * Used with the merge vertices option.
   * \{ */
  
  /* Callbacks for `BM_vert_pair_shared_face_cb` */
  
  struct EDBMSplitBestFaceData {
    BMEdge **edgenet;
    int edgenet_len;
  
    /**
     * Track the range of vertices in edgenet along the faces normal,
     * find the lowest since it's most likely to be most co-planar with the face.
     */
    float best_face_range_on_normal_axis;
    BMFace *r_best_face;
  };
  
  static bool bm_vert_pair_share_best_splittable_face_cb(BMFace *f,
                                                         BMLoop *l_a,
                                                         BMLoop *l_b,
                                                         void *userdata)
  {
    struct EDBMSplitBestFaceData *data = userdata;
    float no[3];
    copy_v3_v3(no, f->no);
  
    float min = dot_v3v3(l_a->v->co, no);
    float max = dot_v3v3(l_b->v->co, no);
    if (min > max) {
      SWAP(float, min, max);
    }
  
    BMVert *v_test = l_b->v;
    BMEdge **e_iter = &data->edgenet[0];
    int verts_len = data->edgenet_len - 1;
    for (int i = verts_len; i--; e_iter++) {
      v_test = BM_edge_other_vert(*e_iter, v_test);
      if (!BM_face_point_inside_test(f, v_test->co)) {
        return false;
      }
      float dot = dot_v3v3(v_test->co, no);
      if (dot < min) {
        min = dot;
      }
      if (dot > max) {
        max = dot;
      }
    }
  
    const float test_face_range_on_normal_axis = max - min;
    if (test_face_range_on_normal_axis < data->best_face_range_on_normal_axis) {
      data->best_face_range_on_normal_axis = test_face_range_on_normal_axis;
      data->r_best_face = f;
    }
  
    return false;
  }
  
  /* find the best splittable face between the two vertices. */
  static bool bm_vert_pair_share_splittable_face_cb(BMFace *UNUSED(f),
                                                    BMLoop *l_a,
                                                    BMLoop *l_b,
                                                    void *userdata)
  {
    float(*data)[3] = userdata;
    float *v_a_co = data[0];
    float *v_a_b_dir = data[1];
  
    float lambda;
    if (isect_ray_seg_v3(v_a_co, v_a_b_dir, l_a->prev->v->co, l_a->next->v->co, &lambda)) {
      if (IN_RANGE(lambda, 0.0f, 1.0f)) {
        return true;
      }
      else if (isect_ray_seg_v3(v_a_co, v_a_b_dir, l_b->prev->v->co, l_b->next->v->co, &lambda)) {
        return IN_RANGE(lambda, 0.0f, 1.0f);
      }
    }
    return false;
  }
  
  void BM_vert_weld_linked_wire_edges_into_linked_faces(
      BMesh *bm, BMVert *v, const float epsilon, BMEdge **r_edgenet[], int *r_edgenet_alloc_len)
  {
    BMEdge **edgenet = *r_edgenet;
    int edgenet_alloc_len = *r_edgenet_alloc_len;
  
    BMIter iter;
    BMEdge *e;
    BM_ITER_ELEM (e, &iter, v, BM_EDGES_OF_VERT) {
      int edgenet_len = 0;
      BMVert *v_other = v;
      while (BM_edge_is_wire(e)) {
        if (edgenet_alloc_len == edgenet_len) {
          edgenet_alloc_len = (edgenet_alloc_len + 1) * 2;
          edgenet = MEM_reallocN(edgenet, (edgenet_alloc_len) * sizeof(*edgenet));
        }
        edgenet[edgenet_len++] = e;
        v_other = BM_edge_other_vert(e, v_other);
        if (v_other == v) {
          /* Endless loop. */
          break;
        }
  
        BMEdge *e_next = BM_DISK_EDGE_NEXT(e, v_other);
        if (e_next == e) {
          /* Vert is wire_endpoint. */
          edgenet_len = 0;
          break;
        }
  
        BMEdge *e_test = e_next;
        while ((e_test = BM_DISK_EDGE_NEXT(e_test, v_other)) != e) {
          if (e_test->l) {
            /* Vert is linked to a face. */
            goto l_break;
          }
        }
  
        e = e_next;
      }
  
      BMLoop *dummy;
      BMFace *best_face;
  
    l_break:
      if (edgenet_len == 0) {
        /* Nothing to do. */
        continue;
      }
      if (edgenet_len == 1) {
        float data[2][3];
        copy_v3_v3(data[0], v_other->co);
        sub_v3_v3v3(data[1], v->co, data[0]);
        best_face = BM_vert_pair_shared_face_cb(
            v_other, v, true, bm_vert_pair_share_splittable_face_cb, &data, &dummy, &dummy);
      }
      else {
        struct EDBMSplitBestFaceData data = {
            .edgenet = edgenet,
            .edgenet_len = edgenet_len,
            .best_face_range_on_normal_axis = FLT_MAX,
            .r_best_face = NULL,
        };
        BM_vert_pair_shared_face_cb(
            v_other, v, true, bm_vert_pair_share_best_splittable_face_cb, &data, &dummy, &dummy);
  
        if (data.r_best_face) {
          float no[3], min = FLT_MAX, max = -FLT_MAX;
          copy_v3_v3(no, data.r_best_face->no);
          BMVert *v_test;
          BMIter f_iter;
          BM_ITER_ELEM (v_test, &f_iter, data.r_best_face, BM_VERTS_OF_FACE) {
            float dot = dot_v3v3(v_test->co, no);
            if (dot < min) {
              min = dot;
            }
            if (dot > max) {
              max = dot;
            }
          }
          float range = max - min + 2 * epsilon;
          if (range < data.best_face_range_on_normal_axis) {
            data.r_best_face = NULL;
          }
        }
        best_face = data.r_best_face;
      }
  
      if (best_face) {
        BM_face_split_edgenet(bm, best_face, edgenet, edgenet_len, NULL, NULL);
      }
    }
  
    *r_edgenet = edgenet;
    *r_edgenet_alloc_len = edgenet_alloc_len;
  }
  
  /** \} */
  
  /* -------------------------------------------------------------------- */
  /** \name Auto-Merge & Split Selection
   *
   * Used after transform operations.
   * \{ */
  
  struct EDBMSplitElem {
    union {
      BMElem *elem;
      BMVert *vert;
      struct {
        BMEdge *edge;
        float lambda;
      };
    };
  };
  
  /* -------------------------------------------------------------------- */
  /* Overlap Callbacks */
  
  struct EDBMSplitData {
    BMesh *bm;
    BLI_Stack *pair_stack;
    int cut_edges_a_len;
    int cut_edges_b_len;
    float dist_sq;
    float dist_sq_sq;
  };
  
  /* Utils */
  
  static void bm_vert_pair_elem_setup_ex(BMVert *v,
                                         float edge_index,
                                         struct EDBMSplitElem *r_pair_elem)
  {
    BLI_assert(v->head.index == -1);
    v->head.index = edge_index;
    r_pair_elem->vert = v;
  }
  
  static void bm_edge_pair_elem_setup(BMEdge *e,
                                      float lambda,
                                      int *r_data_cut_edges_len,
                                      struct EDBMSplitElem *r_pair_elem)
  {
    r_pair_elem->edge = e;
    r_pair_elem->lambda = lambda;
  
    e->head.index++;
    /* Obs: Check Multithread. */
    if (BM_elem_flag_test(e, BM_ELEM_TAG)) {
      BM_elem_flag_disable(e, BM_ELEM_TAG);
      (*r_data_cut_edges_len)++;
    }
  }
  
  /* Util for Vert x Edge and Edge x Edge callbacks */
  static bool bm_vertxedge_isect_impl_ex(BMVert *v,
                                         BMEdge *e,
                                         int edge_index,
                                         const float co[3],
                                         const float dir[3],
                                         float lambda,
                                         float data_dist_sq,
                                         int *data_cut_edges_len,
                                         struct EDBMSplitElem r_pair[2])
  {
    BLI_assert(v->head.index == -1);
  
    BMVert *e_v;
    float dist_sq_vert_factor;
  
    if (lambda < 0.5f) {
      e_v = e->v1;
      dist_sq_vert_factor = lambda;
    }
    else {
      e_v = e->v2;
      dist_sq_vert_factor = 1.0f - lambda;
    }
  
    if (v != e_v) {
      CLAMP(lambda, 0.0f, 1.0f);
  
      float near[3];
      madd_v3_v3v3fl(near, co, dir, lambda);
  
      float dist_sq = len_squared_v3v3(v->co, near);
      if (dist_sq < data_dist_sq) {
        float dist_sq_vert = SQUARE(dist_sq_vert_factor) * len_squared_v3(dir);
        if (dist_sq_vert < data_dist_sq) {
          if (e_v->head.index != -1) {
            /* Vertex already has an intersection. */
            return false;
          }
  
          bm_vert_pair_elem_setup_ex(e_v, -2, &r_pair[1]);
        }
        else {
          bm_edge_pair_elem_setup(e, lambda, data_cut_edges_len, &r_pair[1]);
        }
  
        bm_vert_pair_elem_setup_ex(v, edge_index, &r_pair[0]);
        return true;
      }
    }
  
    return false;
  }
  
  /* Vertex x Vertex Callback */
  
  static bool bm_vertxvert_isect_cb(void *userdata, int index_a, int index_b, int UNUSED(thread))
  {
    struct EDBMSplitData *data = userdata;
    BMVert *v_a = BM_vert_at_index(data->bm, index_a);
    BMVert *v_b = BM_vert_at_index(data->bm, index_b);
  
    struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
  
    BLI_assert(v_a->head.index == -1);
  
    /* Set index -2 for sure that it will not repeat keys in `targetmap`. */
    bm_vert_pair_elem_setup_ex(v_a, -2, &pair[0]);
    bm_vert_pair_elem_setup_ex(v_b, -1, &pair[1]);
  
    return true;
  }
  
  /* Vertex x Edge and Edge x Vertex Callbacks */
  
  static int bm_vertxedge_isect_impl(BMesh *bm,
                                     int vert_index,
                                     int edge_index,
                                     float data_dist_sq,
                                     int *data_cut_edges_len,
                                     struct EDBMSplitElem r_pair[2])
  {
    BMVert *v = BM_vert_at_index(bm, vert_index);
    BMEdge *e = BM_edge_at_index(bm, edge_index);
  
    if (v->head.index != -1) {
      /* Only one vertex per edge. */
      return false;
    }
  
    float co[3], dir[3], lambda;
    copy_v3_v3(co, e->v1->co);
    sub_v3_v3v3(dir, e->v2->co, co);
    lambda = ray_point_factor_v3_ex(v->co, co, dir, 0.0f, -1.0f);
  
    return bm_vertxedge_isect_impl_ex(
        v, e, edge_index, co, dir, lambda, data_dist_sq, data_cut_edges_len, r_pair);
  }
  
  static bool bm_vertxedge_isect_cb(void *userdata, int index_a, int index_b, int UNUSED(thread))
  {
    struct EDBMSplitData *data = userdata;
    struct EDBMSplitElem pair_tmp[2];
    if (bm_vertxedge_isect_impl(
            data->bm, index_a, index_b, data->dist_sq, &data->cut_edges_b_len, pair_tmp)) {
      struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
      pair[0] = pair_tmp[0];
      pair[1] = pair_tmp[1];
  
      return true;
    }
  
    return false;
  }
  
  static bool bm_edgexvert_isect_cb(void *userdata, int index_a, int index_b, int UNUSED(thread))
  {
    struct EDBMSplitData *data = userdata;
    struct EDBMSplitElem pair_tmp[2];
    if (bm_vertxedge_isect_impl(
            data->bm, index_b, index_a, data->dist_sq, &data->cut_edges_a_len, pair_tmp)) {
      struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
      pair[0] = pair_tmp[1];
      pair[1] = pair_tmp[0];
  
      return true;
    }
  
    return false;
  }
  
  /* Edge x Edge Callbacks */
  
  static void bm_edgexedge_isect_impl(struct EDBMSplitData *data,
                                      int index_a,
                                      int index_b,
                                      BMEdge *e_a,
                                      BMEdge *e_b,
                                      const float co_a[3],
                                      const float dir_a[3],
                                      const float co_b[3],
                                      const float dir_b[3],
                                      float lambda_a,
                                      float lambda_b)
  {
    float dist_sq_va_factor, dist_sq_vb_factor;
    BMVert *e_a_v, *e_b_v;
    if (lambda_a < 0.5f) {
      e_a_v = e_a->v1;
      dist_sq_va_factor = lambda_a;
    }
    else {
      e_a_v = e_a->v2;
      dist_sq_va_factor = 1.0f - lambda_a;
    }
  
    if (lambda_b < 0.5f) {
      e_b_v = e_b->v1;
      dist_sq_vb_factor = lambda_b;
    }
    else {
      e_b_v = e_b->v2;
      dist_sq_vb_factor = 1.0f - lambda_b;
    }
  
    if (e_a_v != e_b_v) {
      CLAMP(lambda_a, 0.0f, 1.0f);
      CLAMP(lambda_b, 0.0f, 1.0f);
  
      float near_a[3], near_b[3];
      madd_v3_v3v3fl(near_a, co_a, dir_a, lambda_a);
      madd_v3_v3v3fl(near_b, co_b, dir_b, lambda_b);
  
      float dist_sq = len_squared_v3v3(near_a, near_b);
      if (dist_sq < data->dist_sq) {
        struct EDBMSplitElem pair_tmp[2];
  
        float dist_sq_va = SQUARE(dist_sq_va_factor) * len_squared_v3(dir_a);
        float dist_sq_vb = SQUARE(dist_sq_vb_factor) * len_squared_v3(dir_b);
  
        if (dist_sq_va < data->dist_sq) {
          if (e_a_v->head.index != -1) {
            /* Only one vertex per edge. */
            return;
          }
          bm_vert_pair_elem_setup_ex(e_a_v, index_b, &pair_tmp[0]);
        }
  
        if (dist_sq_vb < data->dist_sq) {
          if (e_b_v->head.index != -1) {
            /* Only one vertex per edge. */
            return;
          }
          bm_vert_pair_elem_setup_ex(e_b_v, index_a, &pair_tmp[1]);
        }
        else {
          bm_edge_pair_elem_setup(e_b, lambda_b, &data->cut_edges_b_len, &pair_tmp[1]);
        }
  
        /* Don't setup edges before a return. */
        if (dist_sq_va >= data->dist_sq) {
          bm_edge_pair_elem_setup(e_a, lambda_a, &data->cut_edges_a_len, &pair_tmp[0]);
        }
  
        struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
        pair[0] = pair_tmp[0];
        pair[1] = pair_tmp[1];
      }
    }
  }
  
  static bool bm_edgexedge_isect_cb(void *userdata, int index_a, int index_b, int UNUSED(thread))
  {
    bool ret = false;
    struct EDBMSplitData *data = userdata;
    BMEdge *e_a = BM_edge_at_index_find(data->bm, index_a);
    BMEdge *e_b = BM_edge_at_index_find(data->bm, index_b);
  
    float co_a[3], dir_a[3], co_b[3], dir_b[3];
    copy_v3_v3(co_a, e_a->v1->co);
    sub_v3_v3v3(dir_a, e_a->v2->co, co_a);
  
    copy_v3_v3(co_b, e_b->v1->co);
    sub_v3_v3v3(dir_b, e_b->v2->co, co_b);
  
    float lambda_a, lambda_b;
    /* Using with dist^4 as `epsilon` is not the best solution, but it fits in most cases. */
    if (isect_ray_ray_epsilon_v3(co_a, dir_a, co_b, dir_b, data->dist_sq_sq, &lambda_a, &lambda_b)) {
      if (ELEM(index_b, e_a->v1->head.index, e_a->v2->head.index) ||
          ELEM(index_a, e_b->v1->head.index, e_b->v2->head.index)) {
        return ret;
      }
  
      /* Edge x Edge returns always false. */
      bm_edgexedge_isect_impl(
          data, index_a, index_b, e_a, e_b, co_a, dir_a, co_b, dir_b, lambda_a, lambda_b);
    }
    else {
      /* Parallel */
      struct EDBMSplitElem pair_tmp[2];
      float vec[3], len_sq_a, len_sq_b, lambda;
      sub_v3_v3v3(vec, co_b, co_a);
      len_sq_a = len_squared_v3(dir_a);
      len_sq_b = len_squared_v3(dir_b);
  
      if (!ELEM(e_b->v1, e_a->v1, e_a->v2) && e_b->v1->head.index == -1) {
        lambda = dot_v3v3(vec, dir_a) / len_sq_a;
        if (bm_vertxedge_isect_impl_ex(e_b->v1,
                                       e_a,
                                       index_a,
                                       co_a,
                                       dir_a,
                                       lambda,
                                       data->dist_sq,
                                       &data->cut_edges_a_len,
                                       pair_tmp)) {
          struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
          pair[0] = pair_tmp[1];
          pair[1] = pair_tmp[0];
          ret |= true;
        }
      }
  
      if (!ELEM(e_a->v1, e_b->v1, e_b->v2) && e_a->v1->head.index == -1) {
        lambda = -dot_v3v3(vec, dir_b) / len_sq_b;
        if (bm_vertxedge_isect_impl_ex(e_a->v1,
                                       e_b,
                                       index_b,
                                       co_b,
                                       dir_b,
                                       lambda,
                                       data->dist_sq,
                                       &data->cut_edges_b_len,
                                       pair_tmp)) {
          struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
          pair[0] = pair_tmp[0];
          pair[1] = pair_tmp[1];
          ret |= true;
        }
      }
  
      add_v3_v3(vec, dir_b);
      if (!ELEM(e_b->v2, e_a->v1, e_a->v2) && e_b->v2->head.index == -1) {
        lambda = dot_v3v3(vec, dir_a) / len_sq_a;
        if (bm_vertxedge_isect_impl_ex(e_b->v2,
                                       e_a,
                                       index_a,
                                       co_a,
                                       dir_a,
                                       lambda,
                                       data->dist_sq,
                                       &data->cut_edges_a_len,
                                       pair_tmp)) {
          struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
          pair[0] = pair_tmp[1];
          pair[1] = pair_tmp[0];
          ret |= true;
        }
      }
  
      sub_v3_v3(vec, dir_a);
      if (!ELEM(e_a->v2, e_b->v1, e_b->v2) && e_a->v2->head.index == -1) {
        lambda = 1.0f - dot_v3v3(vec, dir_b) / len_sq_b;
        if (bm_vertxedge_isect_impl_ex(e_a->v2,
                                       e_b,
                                       index_b,
                                       co_b,
                                       dir_b,
                                       lambda,
                                       data->dist_sq,
                                       &data->cut_edges_b_len,
                                       pair_tmp)) {
          struct EDBMSplitElem *pair = BLI_stack_push_r(data->pair_stack);
          pair[0] = pair_tmp[0];
          pair[1] = pair_tmp[1];
          ret |= true;
        }
      }
    }
  
    return ret;
  }
  
  /* -------------------------------------------------------------------- */
  /* BVHTree Overlap Function */
  
  static void bvhtree_overlap_thread_safe(const BVHTree *tree1,
                                          const BVHTree *tree2,
                                          BVHTree_OverlapCallback callback,
                                          void *userdata)
  {
    BLI_bvhtree_overlap_ex(tree1, tree2, NULL, callback, userdata, 1, 0);
  }
  
  /* -------------------------------------------------------------------- */
  /* Callbacks for `BLI_qsort_r` */
  
  static int sort_cmp_by_lambda_a_cb(const void *index1_v, const void *index2_v, void *keys_v)
  {
    const struct EDBMSplitElem(*pair_array)[2] = keys_v;
    const int index1 = *(int *)index1_v;
    const int index2 = *(int *)index2_v;
  
    if (pair_array[index1][0].lambda > pair_array[index2][0].lambda) {
      return 1;
    }
    else {
      return -1;
    }
  }
  
  static int sort_cmp_by_lambda_b_cb(const void *index1_v, const void *index2_v, void *keys_v)
  {
    const struct EDBMSplitElem(*pair_array)[2] = keys_v;
    const int index1 = *(int *)index1_v;
    const int index2 = *(int *)index2_v;
  
    if (pair_array[index1][1].lambda > pair_array[index2][1].lambda) {
      return 1;
    }
    else {
      return -1;
    }
  }
  
  /* -------------------------------------------------------------------- */
  /* Main API */
  
  bool BM_mesh_intersect_edges(BMesh *bm, const char hflag, const float dist, GHash *r_targetmap)
  {
    bool ok = false;
  
    BMIter iter;
    BMVert *v;
    BMEdge *e;
    int i;
  
    BM_mesh_elem_table_ensure(bm, BM_VERT | BM_EDGE);
  
    /* Store all intersections in this array. */
    struct EDBMSplitElem(*pair_iter)[2], (*pair_array)[2] = NULL;
    BLI_Stack *pair_stack = BLI_stack_new(sizeof(*pair_array), __func__);
    int pair_len = 0;
  
    float dist_sq = SQUARE(dist);
    struct EDBMSplitData data = {
        .bm = bm,
        .pair_stack = pair_stack,
        .cut_edges_a_len = 0,
        .cut_edges_b_len = 0,
        .dist_sq = dist_sq,
        .dist_sq_sq = SQUARE(dist_sq),
    };
  
    /* tag and count the verts to be tested. */
    int verts_act_len = 0, verts_remain_len = 0;
    int loose_verts_act_len = 0, loose_verts_remain_len = 0;
    BM_ITER_MESH (v, &iter, bm, BM_VERTS_OF_MESH) {
      if (BM_elem_flag_test(v, hflag)) {
        BM_elem_flag_enable(v, BM_ELEM_TAG);
        v->head.index = -1;
        verts_act_len++;
        if (!v->e) {
          loose_verts_act_len++;
        }
      }
      else {
        BM_elem_flag_disable(v, BM_ELEM_TAG);
        if (!BM_elem_flag_test(v, BM_ELEM_HIDDEN)) {
          v->head.index = -1;
          verts_remain_len++;
          if (!v->e) {
            loose_verts_remain_len++;
          }
        }
      }
    }
    bm->elem_index_dirty |= BM_VERT;
  
    /* Start the creation of BVHTrees. */
    BVHTree *tree_loose_verts_act = NULL, *tree_loose_verts_remain = NULL;
    if (loose_verts_act_len) {
      tree_loose_verts_act = BLI_bvhtree_new(loose_verts_act_len, dist, 2, KDOP_AXIS_LEN);
    }
    if (loose_verts_remain_len) {
      tree_loose_verts_remain = BLI_bvhtree_new(loose_verts_remain_len, 0.0f, 2, KDOP_AXIS_LEN);
    }
  
    if (tree_loose_verts_act || tree_loose_verts_remain) {
      BM_ITER_MESH_INDEX (v, &iter, bm, BM_VERTS_OF_MESH, i) {
        if (BM_elem_flag_test(v, BM_ELEM_TAG)) {
          if (tree_loose_verts_act && !v->e) {
            BLI_bvhtree_insert(tree_loose_verts_act, i, v->co, 1);
          }
        }
        else if (tree_loose_verts_remain && !v->e && !BM_elem_flag_test(v, BM_ELEM_HIDDEN)) {
          BLI_bvhtree_insert(tree_loose_verts_remain, i, v->co, 1);
        }
      }
      if (tree_loose_verts_act) {
        BLI_bvhtree_balance(tree_loose_verts_act);
      }
  
      if (tree_loose_verts_remain) {
        BLI_bvhtree_balance(tree_loose_verts_remain);
      }
  
      if (tree_loose_verts_act && tree_loose_verts_remain) {
        /* First pair search. */
        bvhtree_overlap_thread_safe(
            tree_loose_verts_act, tree_loose_verts_remain, bm_vertxvert_isect_cb, &data);
      }
    }
  
    /* Tag and count the edges. */
    int edges_act_len = 0, edges_remain_len = 0;
    BM_ITER_MESH (e, &iter, bm, BM_EDGES_OF_MESH) {
      if (BM_elem_flag_test(e->v1, BM_ELEM_TAG) || BM_elem_flag_test(e->v2, BM_ELEM_TAG)) {
        BM_elem_flag_enable(e, BM_ELEM_TAG);
        edges_act_len++;
      }
      else {
        BM_elem_flag_disable(e, BM_ELEM_TAG);
        if (!BM_elem_flag_test(e, BM_ELEM_HIDDEN)) {
          edges_remain_len++;
        }
      }
    }
  
    if (edges_remain_len) {
      BVHTree *tree_edges_act = NULL, *tree_edges_remain = NULL;
      tree_edges_remain = BLI_bvhtree_new(edges_remain_len, 0.0f, 2, KDOP_AXIS_LEN);
      if (edges_act_len) {
        tree_edges_act = BLI_bvhtree_new(edges_act_len, dist, 2, KDOP_AXIS_LEN);
      }
  
      BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
        float co[2][3];
        if (BM_elem_flag_test(e, BM_ELEM_TAG)) {
          if (tree_edges_act) {
            e->head.index = 0;
            copy_v3_v3(co[0], e->v1->co);
            copy_v3_v3(co[1], e->v2->co);
            BLI_bvhtree_insert(tree_edges_act, i, co[0], 2);
          }
        }
        else if (!BM_elem_flag_test(e, BM_ELEM_HIDDEN)) {
          /* Tag used in the overlap callbacks. */
          BM_elem_flag_enable(e, BM_ELEM_TAG);
          e->head.index = 0;
          copy_v3_v3(co[0], e->v1->co);
          copy_v3_v3(co[1], e->v2->co);
          BLI_bvhtree_insert(tree_edges_remain, i, co[0], 2);
        }
      }
      /* Use `e->head.index` to count intersections. */
      bm->elem_index_dirty |= BM_EDGE;
  
      BLI_bvhtree_balance(tree_edges_remain);
      if (tree_edges_act) {
        BLI_bvhtree_balance(tree_edges_act);
      }
  
      if (tree_edges_act) {
        /* Edge x Edge */
        bvhtree_overlap_thread_safe(tree_edges_act, tree_edges_remain, bm_edgexedge_isect_cb, &data);
  
        if (tree_loose_verts_remain) {
          /* Edge x Vert */
          bvhtree_overlap_thread_safe(
              tree_edges_act, tree_loose_verts_remain, bm_edgexvert_isect_cb, &data);
        }
  
        BLI_bvhtree_free(tree_edges_act);
      }
  
      if (tree_loose_verts_act) {
        /* Vert x Edge */
        bvhtree_overlap_thread_safe(
            tree_loose_verts_act, tree_edges_remain, bm_vertxedge_isect_cb, &data);
      }
  
      BLI_bvhtree_free(tree_edges_remain);
  
      pair_len = BLI_stack_count(pair_stack);
      if (pair_len) {
        pair_array = MEM_mallocN(sizeof(*pair_array) * pair_len, __func__);
        BLI_stack_pop_n_reverse(pair_stack, pair_array, pair_len);
  
        /* Map intersections per edge. */
        union {
          struct {
            int cuts_len;
            int cuts_index[];
          };
          int as_int[0];
        } * e_map_iter, *e_map;
  
        size_t e_map_size = (max_ii(data.cut_edges_a_len, data.cut_edges_b_len) * sizeof(*e_map)) +
                            (pair_len * sizeof(*(e_map->cuts_index)));
  
        e_map = MEM_mallocN(e_map_size, __func__);
  
        /* Convert every pair to Vert x Vert. */
        for (int pair = 0; pair < 2; pair++) {
          int map_len = 0;
          pair_iter = &pair_array[0];
          for (i = 0; i < pair_len; i++, pair_iter++) {
            if ((*pair_iter)[pair].elem->head.htype != BM_EDGE) {
              /* Take the opportunity to set all vert indices to -1 again. */
              (*pair_iter)[pair].elem->head.index = -1;
              continue;
            }
            e = (*pair_iter)[pair].edge;
            if (!BM_elem_flag_test(e, BM_ELEM_TAG)) {
              BM_elem_flag_enable(e, BM_ELEM_TAG);
              int e_cuts_len = e->head.index;
  
              e_map_iter = (void *)&e_map->as_int[map_len];
              e_map_iter->cuts_len = e_cuts_len;
              e_map_iter->cuts_index[0] = i;
  
              /* Use `e->head.index` to indicate which slot to fill with the `cut` index. */
              e->head.index = map_len + 1;
              map_len += 1 + e_cuts_len;
            }
            else {
              e_map->as_int[++e->head.index] = i;
            }
          }
  
          /* Split Edges A to set all Vert x Edge. */
          for (i = 0; i < map_len;
               e_map_iter = (void *)&e_map->as_int[i], i += 1 + e_map_iter->cuts_len) {
  
            /* sort by lambda. */
            BLI_qsort_r(e_map_iter->cuts_index,
                        e_map_iter->cuts_len,
                        sizeof(*(e_map->cuts_index)),
                        pair == 0 ? sort_cmp_by_lambda_a_cb : sort_cmp_by_lambda_b_cb,
                        pair_array);
  
            float lambda, lambda_prev = 0.0f;
            for (int j = 0; j < e_map_iter->cuts_len; j++) {
              struct EDBMSplitElem *pair_elem = &pair_array[e_map_iter->cuts_index[j]][pair];
              lambda = (pair_elem->lambda - lambda_prev) / (1.0f - lambda_prev);
              lambda_prev = pair_elem->lambda;
              e = pair_elem->edge;
  
              BMVert *v_new = BM_edge_split(bm, e, e->v1, NULL, lambda);
              v_new->head.index = -1;
              pair_elem->vert = v_new;
            }
          }
        }
  
        MEM_freeN(e_map);
      }
    }
  
    BLI_bvhtree_free(tree_loose_verts_act);
    BLI_bvhtree_free(tree_loose_verts_remain);
  
    if (r_targetmap) {
      if (pair_array == NULL) {
        pair_len = BLI_stack_count(pair_stack);
        if (pair_len) {
          pair_array = MEM_mallocN(sizeof(*pair_array) * pair_len, __func__);
          BLI_stack_pop_n_reverse(pair_stack, pair_array, pair_len);
        }
      }
  
      if (pair_array) {
        /* Organize the vertices in the order they will be merged. */
        pair_iter = &pair_array[0];
        for (i = 0; i < pair_len; i++, pair_iter++) {
          BLI_assert((*pair_iter)[0].elem->head.htype == BM_VERT);
          BLI_assert((*pair_iter)[1].elem->head.htype == BM_VERT);
          BLI_assert((*pair_iter)[0].elem != (*pair_iter)[1].elem);
  
          BLI_ghash_insert(r_targetmap, (*pair_iter)[0].vert, (*pair_iter)[1].vert);
        }
  
        ok = true;
      }
    }
  
    BLI_stack_free(pair_stack);
    if (pair_array) {
      MEM_freeN(pair_array);
    }
  
    return ok;
  }
  
  /** \} */
  