test2/source/blender/blenkernel/intern/mesh_remesh_voxel.cc

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

/** \file
 * \ingroup bke
 */

#include <cctype>
#include <cfloat>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctime>

#include "MEM_guardedalloc.h"

#include "BLI_array.hh"
#include "BLI_array_utils.hh"
#include "BLI_enumerable_thread_specific.hh"
#include "BLI_index_range.hh"
#include "BLI_math_vector.h"
#include "BLI_span.hh"
#include "BLI_task.hh"

#include "DNA_mesh_types.h"
#include "DNA_meshdata_types.h"

#include "BKE_attribute.h"
#include "BKE_attribute.hh"
#include "BKE_attribute_math.hh"
#include "BKE_bvhutils.hh"
#include "BKE_customdata.hh"
#include "BKE_editmesh.hh"
#include "BKE_lib_id.h"
#include "BKE_mesh.hh"
#include "BKE_mesh_mapping.hh"
#include "BKE_mesh_remesh_voxel.hh" /* own include */
#include "BKE_mesh_runtime.hh"
#include "BKE_mesh_sample.hh"

#include "bmesh_tools.h"

#ifdef WITH_OPENVDB
#  include <openvdb/openvdb.h>
#  include <openvdb/tools/MeshToVolume.h>
#  include <openvdb/tools/VolumeToMesh.h>
#endif

#ifdef WITH_QUADRIFLOW
#  include "quadriflow_capi.hpp"
#endif

using blender::Array;
using blender::float3;
using blender::IndexRange;
using blender::MutableSpan;
using blender::Span;

#ifdef WITH_QUADRIFLOW
static Mesh *remesh_quadriflow(const Mesh *input_mesh,
                               int target_faces,
                               int seed,
                               bool preserve_sharp,
                               bool preserve_boundary,
                               bool adaptive_scale,
                               void (*update_cb)(void *, float progress, int *cancel),
                               void *update_cb_data)
{
  const Span<float3> input_positions = input_mesh->vert_positions();
  const Span<int> input_corner_verts = input_mesh->corner_verts();
  const Span<MLoopTri> looptris = input_mesh->looptris();

  /* Gather the required data for export to the internal quadriflow mesh format. */
  Array<MVertTri> verttri(looptris.size());
  BKE_mesh_runtime_verttri_from_looptri(
      verttri.data(), input_corner_verts.data(), looptris.data(), looptris.size());

  const int totfaces = looptris.size();
  const int totverts = input_mesh->totvert;
  Array<int> faces(totfaces * 3);

  for (const int i : IndexRange(totfaces)) {
    MVertTri &vt = verttri[i];
    faces[i * 3] = vt.tri[0];
    faces[i * 3 + 1] = vt.tri[1];
    faces[i * 3 + 2] = vt.tri[2];
  }

  /* Fill out the required input data */
  QuadriflowRemeshData qrd;

  qrd.totfaces = totfaces;
  qrd.totverts = totverts;
  qrd.verts = (float *)input_positions.data();
  qrd.faces = faces.data();
  qrd.target_faces = target_faces;

  qrd.preserve_sharp = preserve_sharp;
  qrd.preserve_boundary = preserve_boundary;
  qrd.adaptive_scale = adaptive_scale;
  qrd.minimum_cost_flow = false;
  qrd.aggresive_sat = false;
  qrd.rng_seed = seed;

  qrd.out_faces = nullptr;

  /* Run the remesher */
  QFLOW_quadriflow_remesh(&qrd, update_cb, update_cb_data);

  if (qrd.out_faces == nullptr) {
    /* The remeshing was canceled */
    return nullptr;
  }

  if (qrd.out_totfaces == 0) {
    /* Meshing failed */
    MEM_freeN(qrd.out_faces);
    MEM_freeN(qrd.out_verts);
    return nullptr;
  }

  /* Construct the new output mesh */
  Mesh *mesh = BKE_mesh_new_nomain(qrd.out_totverts, 0, qrd.out_totfaces, qrd.out_totfaces * 4);
  BKE_mesh_copy_parameters(mesh, input_mesh);
  MutableSpan<int> face_offsets = mesh->face_offsets_for_write();
  MutableSpan<int> corner_verts = mesh->corner_verts_for_write();

  blender::offset_indices::fill_constant_group_size(4, 0, face_offsets);

  mesh->vert_positions_for_write().copy_from(
      Span(reinterpret_cast<float3 *>(qrd.out_verts), qrd.out_totverts));

  for (const int i : IndexRange(qrd.out_totfaces)) {
    const int loopstart = i * 4;
    corner_verts[loopstart] = qrd.out_faces[loopstart];
    corner_verts[loopstart + 1] = qrd.out_faces[loopstart + 1];
    corner_verts[loopstart + 2] = qrd.out_faces[loopstart + 2];
    corner_verts[loopstart + 3] = qrd.out_faces[loopstart + 3];
  }

  BKE_mesh_calc_edges(mesh, false, false);

  MEM_freeN(qrd.out_faces);
  MEM_freeN(qrd.out_verts);

  return mesh;
}
#endif

Mesh *BKE_mesh_remesh_quadriflow(const Mesh *mesh,
                                 int target_faces,
                                 int seed,
                                 bool preserve_sharp,
                                 bool preserve_boundary,
                                 bool adaptive_scale,
                                 void (*update_cb)(void *, float progress, int *cancel),
                                 void *update_cb_data)
{
#ifdef WITH_QUADRIFLOW
  if (target_faces <= 0) {
    target_faces = -1;
  }
  return remesh_quadriflow(mesh,
                           target_faces,
                           seed,
                           preserve_sharp,
                           preserve_boundary,
                           adaptive_scale,
                           update_cb,
                           update_cb_data);
#else
  UNUSED_VARS(mesh,
              target_faces,
              seed,
              preserve_sharp,
              preserve_boundary,
              adaptive_scale,
              update_cb,
              update_cb_data);
  return nullptr;
#endif
}

#ifdef WITH_OPENVDB
static openvdb::FloatGrid::Ptr remesh_voxel_level_set_create(const Mesh *mesh,
                                                             const float voxel_size)
{
  const Span<float3> positions = mesh->vert_positions();
  const Span<int> corner_verts = mesh->corner_verts();
  const Span<MLoopTri> looptris = mesh->looptris();

  std::vector<openvdb::Vec3s> points(mesh->totvert);
  std::vector<openvdb::Vec3I> triangles(looptris.size());

  for (const int i : IndexRange(mesh->totvert)) {
    const float3 &co = positions[i];
    points[i] = openvdb::Vec3s(co.x, co.y, co.z);
  }

  for (const int i : IndexRange(looptris.size())) {
    const MLoopTri &loop_tri = looptris[i];
    triangles[i] = openvdb::Vec3I(corner_verts[loop_tri.tri[0]],
                                  corner_verts[loop_tri.tri[1]],
                                  corner_verts[loop_tri.tri[2]]);
  }

  openvdb::math::Transform::Ptr transform = openvdb::math::Transform::createLinearTransform(
      voxel_size);
  openvdb::FloatGrid::Ptr grid = openvdb::tools::meshToLevelSet<openvdb::FloatGrid>(
      *transform, points, triangles, 1.0f);

  return grid;
}

static Mesh *remesh_voxel_volume_to_mesh(const openvdb::FloatGrid::Ptr level_set_grid,
                                         const float isovalue,
                                         const float adaptivity,
                                         const bool relax_disoriented_triangles)
{
  std::vector<openvdb::Vec3s> vertices;
  std::vector<openvdb::Vec4I> quads;
  std::vector<openvdb::Vec3I> tris;
  openvdb::tools::volumeToMesh<openvdb::FloatGrid>(
      *level_set_grid, vertices, tris, quads, isovalue, adaptivity, relax_disoriented_triangles);

  Mesh *mesh = BKE_mesh_new_nomain(
      vertices.size(), 0, quads.size() + tris.size(), quads.size() * 4 + tris.size() * 3);
  MutableSpan<float3> vert_positions = mesh->vert_positions_for_write();
  MutableSpan<int> face_offsets = mesh->face_offsets_for_write();
  MutableSpan<int> mesh_corner_verts = mesh->corner_verts_for_write();

  const int triangle_loop_start = quads.size() * 4;
  if (!face_offsets.is_empty()) {
    blender::offset_indices::fill_constant_group_size(4, 0, face_offsets.take_front(quads.size()));
    blender::offset_indices::fill_constant_group_size(
        3, triangle_loop_start, face_offsets.drop_front(quads.size()));
  }

  for (const int i : vert_positions.index_range()) {
    vert_positions[i] = float3(vertices[i].x(), vertices[i].y(), vertices[i].z());
  }

  for (const int i : IndexRange(quads.size())) {
    const int loopstart = i * 4;
    mesh_corner_verts[loopstart] = quads[i][0];
    mesh_corner_verts[loopstart + 1] = quads[i][3];
    mesh_corner_verts[loopstart + 2] = quads[i][2];
    mesh_corner_verts[loopstart + 3] = quads[i][1];
  }

  for (const int i : IndexRange(tris.size())) {
    const int loopstart = triangle_loop_start + i * 3;
    mesh_corner_verts[loopstart] = tris[i][2];
    mesh_corner_verts[loopstart + 1] = tris[i][1];
    mesh_corner_verts[loopstart + 2] = tris[i][0];
  }

  BKE_mesh_calc_edges(mesh, false, false);

  return mesh;
}
#endif

Mesh *BKE_mesh_remesh_voxel(const Mesh *mesh,
                            const float voxel_size,
                            const float adaptivity,
                            const float isovalue)
{
#ifdef WITH_OPENVDB
  openvdb::FloatGrid::Ptr level_set = remesh_voxel_level_set_create(mesh, voxel_size);
  Mesh *result = remesh_voxel_volume_to_mesh(level_set, isovalue, adaptivity, false);
  BKE_mesh_copy_parameters(result, mesh);
  return result;
#else
  UNUSED_VARS(mesh, voxel_size, adaptivity, isovalue);
  return nullptr;
#endif
}

namespace blender::bke {

static void calc_edge_centers(const Span<float3> positions,
                              const Span<int2> edges,
                              MutableSpan<float3> edge_centers)
{
  for (const int i : edges.index_range()) {
    edge_centers[i] = math::midpoint(positions[edges[i][0]], positions[edges[i][1]]);
  }
}

static void calc_face_centers(const Span<float3> positions,
                              const OffsetIndices<int> faces,
                              const Span<int> corner_verts,
                              MutableSpan<float3> face_centers)
{
  for (const int i : faces.index_range()) {
    face_centers[i] = mesh::face_center_calc(positions, corner_verts.slice(faces[i]));
  }
}

static void find_nearest_tris(const Span<float3> positions,
                              BVHTreeFromMesh &bvhtree,
                              MutableSpan<int> tris)
{
  for (const int i : positions.index_range()) {
    BVHTreeNearest nearest;
    nearest.index = -1;
    nearest.dist_sq = FLT_MAX;
    BLI_bvhtree_find_nearest(
        bvhtree.tree, positions[i], &nearest, bvhtree.nearest_callback, &bvhtree);
    tris[i] = nearest.index;
  }
}

static void find_nearest_tris_parallel(const Span<float3> positions,
                                       BVHTreeFromMesh &bvhtree,
                                       MutableSpan<int> tris)
{
  threading::parallel_for(tris.index_range(), 512, [&](const IndexRange range) {
    find_nearest_tris(positions.slice(range), bvhtree, tris.slice(range));
  });
}

static void find_nearest_verts(const Span<float3> positions,
                               const Span<int> corner_verts,
                               const Span<MLoopTri> src_tris,
                               const Span<float3> dst_positions,
                               const Span<int> nearest_vert_tris,
                               MutableSpan<int> nearest_verts)
{
  threading::parallel_for(dst_positions.index_range(), 512, [&](const IndexRange range) {
    for (const int dst_vert : range) {
      const float3 &dst_position = dst_positions[dst_vert];
      const MLoopTri &src_tri = src_tris[nearest_vert_tris[dst_vert]];

      std::array<float, 3> distances;
      for (const int i : IndexRange(3)) {
        const int src_vert = corner_verts[src_tri.tri[i]];
        distances[i] = math::distance_squared(positions[src_vert], dst_position);
      }

      const int min = std::min_element(distances.begin(), distances.end()) - distances.begin();
      nearest_verts[dst_vert] = corner_verts[src_tri.tri[min]];
    }
  });
}

static void find_nearest_faces(const Span<int> src_tri_faces,
                               const Span<float3> dst_positions,
                               const OffsetIndices<int> dst_faces,
                               const Span<int> dst_corner_verts,
                               BVHTreeFromMesh &bvhtree,
                               MutableSpan<int> nearest_faces)
{
  struct TLS {
    Vector<float3> face_centers;
    Vector<int> tri_indices;
  };
  threading::EnumerableThreadSpecific<TLS> all_tls;
  threading::parallel_for(dst_faces.index_range(), 512, [&](const IndexRange range) {
    TLS &tls = all_tls.local();
    Vector<float3> &face_centers = tls.face_centers;
    face_centers.reinitialize(range.size());
    calc_face_centers(dst_positions, dst_faces.slice(range), dst_corner_verts, face_centers);

    Vector<int> &tri_indices = tls.tri_indices;
    tri_indices.reinitialize(range.size());
    find_nearest_tris(face_centers, bvhtree, tri_indices);

    array_utils::gather(src_tri_faces, tri_indices.as_span(), nearest_faces.slice(range));
  });
}

static void find_nearest_corners(const Span<float3> src_positions,
                                 const OffsetIndices<int> src_faces,
                                 const Span<int> src_corner_verts,
                                 const Span<int> src_tri_faces,
                                 const Span<float3> dst_positions,
                                 const Span<int> dst_corner_verts,
                                 const Span<int> nearest_vert_tris,
                                 MutableSpan<int> nearest_corners)
{
  threading::parallel_for(nearest_corners.index_range(), 512, [&](const IndexRange range) {
    Vector<float, 64> distances;
    for (const int dst_corner : range) {
      const int dst_vert = dst_corner_verts[dst_corner];
      const float3 &dst_position = dst_positions[dst_vert];

      const int src_tri = nearest_vert_tris[dst_vert];
      const IndexRange src_face = src_faces[src_tri_faces[src_tri]];
      const Span<int> src_face_verts = src_corner_verts.slice(src_face);

      /* Find the corner in the face that's closest in the closest face. */
      distances.reinitialize(src_face_verts.size());
      for (const int i : src_face_verts.index_range()) {
        const int src_vert = src_face_verts[i];
        distances[i] = math::distance_squared(src_positions[src_vert], dst_position);
      }

      const int min = std::min_element(distances.begin(), distances.end()) - distances.begin();
      nearest_corners[dst_corner] = src_face[min];
    }
  });
}

static void find_nearest_edges(const Span<float3> src_positions,
                               const Span<int2> src_edges,
                               const OffsetIndices<int> src_faces,
                               const Span<int> src_corner_edges,
                               const Span<int> src_tri_faces,
                               const Span<float3> dst_positions,
                               const Span<int2> dst_edges,
                               BVHTreeFromMesh &bvhtree,
                               MutableSpan<int> nearest_edges)
{
  struct TLS {
    Vector<float3> edge_centers;
    Vector<int> tri_indices;
    Vector<int> face_indices;
    Vector<float> distances;
  };
  threading::EnumerableThreadSpecific<TLS> all_tls;
  threading::parallel_for(nearest_edges.index_range(), 512, [&](const IndexRange range) {
    TLS &tls = all_tls.local();
    Vector<float3> &edge_centers = tls.edge_centers;
    edge_centers.reinitialize(range.size());
    calc_edge_centers(dst_positions, dst_edges.slice(range), edge_centers);

    Vector<int> &tri_indices = tls.tri_indices;
    tri_indices.reinitialize(range.size());
    find_nearest_tris_parallel(edge_centers, bvhtree, tri_indices);

    Vector<int> &face_indices = tls.face_indices;
    face_indices.reinitialize(range.size());
    array_utils::gather(src_tri_faces, tri_indices.as_span(), face_indices.as_mutable_span());

    /* Find the source edge that's closest to the destination edge in the nearest face. Search
     * through the whole face instead of just the triangle because the triangle has edges that
     * might not be actual mesh edges. */
    Vector<float, 64> distances;
    for (const int i : range.index_range()) {
      const int dst_edge = range[i];
      const float3 &dst_position = edge_centers[i];

      const int src_face = face_indices[i];
      const Span<int> src_face_edges = src_corner_edges.slice(src_faces[src_face]);

      distances.reinitialize(src_face_edges.size());
      for (const int i : src_face_edges.index_range()) {
        const int2 src_edge = src_edges[src_face_edges[i]];
        const float3 src_center = math::midpoint(src_positions[src_edge[0]],
                                                 src_positions[src_edge[1]]);
        distances[i] = math::distance_squared(src_center, dst_position);
      }

      const int min = std::min_element(distances.begin(), distances.end()) - distances.begin();
      nearest_edges[dst_edge] = src_face_edges[min];
    }
  });
}

static void gather_attributes(const Span<AttributeIDRef> ids,
                              const AttributeAccessor src_attributes,
                              const eAttrDomain domain,
                              const Span<int> index_map,
                              MutableAttributeAccessor dst_attributes)
{
  for (const AttributeIDRef &id : ids) {
    const GVArraySpan src = *src_attributes.lookup(id, domain);
    const eCustomDataType type = cpp_type_to_custom_data_type(src.type());
    GSpanAttributeWriter dst = dst_attributes.lookup_or_add_for_write_only_span(id, domain, type);
    attribute_math::gather(src, index_map, dst.span);
    dst.finish();
  }
}

void mesh_remesh_reproject_attributes(const Mesh &src, Mesh &dst)
{
  /* Gather attributes to transfer for each domain. This makes it possible to skip
   * building index maps and even the main BVH tree if there are no attributes. */
  const AttributeAccessor src_attributes = src.attributes();
  Vector<AttributeIDRef> point_ids;
  Vector<AttributeIDRef> edge_ids;
  Vector<AttributeIDRef> face_ids;
  Vector<AttributeIDRef> corner_ids;
  src_attributes.for_all([&](const AttributeIDRef &id, const AttributeMetaData &meta_data) {
    if (ELEM(id.name(), "position", ".edge_verts", ".corner_vert", ".corner_edge")) {
      return true;
    }
    switch (meta_data.domain) {
      case ATTR_DOMAIN_POINT:
        point_ids.append(id);
        break;
      case ATTR_DOMAIN_EDGE:
        edge_ids.append(id);
        break;
      case ATTR_DOMAIN_FACE:
        face_ids.append(id);
        break;
      case ATTR_DOMAIN_CORNER:
        corner_ids.append(id);
        break;
      default:
        BLI_assert_unreachable();
        break;
    }
    return true;
  });

  if (point_ids.is_empty() && edge_ids.is_empty() && face_ids.is_empty() && corner_ids.is_empty())
  {
    return;
  }

  const Span<float3> src_positions = src.vert_positions();
  const OffsetIndices src_faces = src.faces();
  const Span<int> src_corner_verts = src.corner_verts();
  const Span<MLoopTri> src_tris = src.looptris();

  /* The main idea in the following code is to trade some complexity in sampling for the benefit of
   * only using and building a single BVH tree. Since sculpt mode doesn't generally deal with loose
   * vertices and edges, we use the standard "triangles" BVH which won't contain them. Also, only
   * relying on a single BVH should reduce memory usage, and work better if the BVH and PBVH are
   * ever merged.
   *
   * One key decision is separating building transfer index maps from actually transferring any
   * attribute data. This is important to keep attribute storage independent from the specifics of
   * the decisions made here, which mainly results in easier refactoring, more generic code, and
   * possibly improved performance from lower cache usage in the "complex" sampling part of the
   * algorithm and the copying itself. */
  BVHTreeFromMesh bvhtree{};
  BKE_bvhtree_from_mesh_get(&bvhtree, &src, BVHTREE_FROM_LOOPTRI, 2);

  const Span<float3> dst_positions = dst.vert_positions();
  const OffsetIndices dst_faces = dst.faces();
  const Span<int> dst_corner_verts = dst.corner_verts();

  MutableAttributeAccessor dst_attributes = dst.attributes_for_write();

  if (!point_ids.is_empty() || !corner_ids.is_empty()) {
    Array<int> vert_nearest_tris(dst_positions.size());
    find_nearest_tris_parallel(dst_positions, bvhtree, vert_nearest_tris);

    if (!point_ids.is_empty()) {
      Array<int> map(dst.totvert);
      find_nearest_verts(
          src_positions, src_corner_verts, src_tris, dst_positions, vert_nearest_tris, map);
      gather_attributes(point_ids, src_attributes, ATTR_DOMAIN_POINT, map, dst_attributes);
    }

    if (!corner_ids.is_empty()) {
      const Span<int> src_tri_faces = src.looptri_faces();
      Array<int> map(dst.totloop);
      find_nearest_corners(src_positions,
                           src_faces,
                           src_corner_verts,
                           src_tri_faces,
                           dst_positions,
                           dst_corner_verts,
                           vert_nearest_tris,
                           map);
      gather_attributes(corner_ids, src_attributes, ATTR_DOMAIN_CORNER, map, dst_attributes);
    }
  }

  if (!edge_ids.is_empty()) {
    const Span<int2> src_edges = src.edges();
    const Span<int> src_corner_edges = src.corner_edges();
    const Span<int> src_tri_faces = src.looptri_faces();
    const Span<int2> dst_edges = dst.edges();
    Array<int> map(dst.totedge);
    find_nearest_edges(src_positions,
                       src_edges,
                       src_faces,
                       src_corner_edges,
                       src_tri_faces,
                       dst_positions,
                       dst_edges,
                       bvhtree,
                       map);
    gather_attributes(edge_ids, src_attributes, ATTR_DOMAIN_EDGE, map, dst_attributes);
  }

  if (!face_ids.is_empty()) {
    const Span<int> src_tri_faces = src.looptri_faces();
    Array<int> map(dst.faces_num);
    find_nearest_faces(src_tri_faces, dst_positions, dst_faces, dst_corner_verts, bvhtree, map);
    gather_attributes(face_ids, src_attributes, ATTR_DOMAIN_FACE, map, dst_attributes);
  }

  if (src.active_color_attribute) {
    BKE_id_attributes_active_color_set(&dst.id, src.active_color_attribute);
  }
  if (src.default_color_attribute) {
    BKE_id_attributes_default_color_set(&dst.id, src.default_color_attribute);
  }

  free_bvhtree_from_mesh(&bvhtree);
}

}  // namespace blender::bke

Mesh *BKE_mesh_remesh_voxel_fix_poles(const Mesh *mesh)
{
  const BMAllocTemplate allocsize = BMALLOC_TEMPLATE_FROM_ME(mesh);

  BMeshCreateParams bmesh_create_params{};
  bmesh_create_params.use_toolflags = true;
  BMesh *bm = BM_mesh_create(&allocsize, &bmesh_create_params);

  BMeshFromMeshParams bmesh_from_mesh_params{};
  bmesh_from_mesh_params.calc_face_normal = true;
  bmesh_from_mesh_params.calc_vert_normal = true;
  BM_mesh_bm_from_me(bm, mesh, &bmesh_from_mesh_params);

  BMVert *v;
  BMEdge *ed, *ed_next;
  BMFace *f, *f_next;
  BMIter iter_a, iter_b;

  /* Merge 3 edge poles vertices that exist in the same face */
  BM_mesh_elem_hflag_disable_all(bm, BM_VERT | BM_EDGE | BM_FACE, BM_ELEM_TAG, false);
  BM_ITER_MESH_MUTABLE (f, f_next, &iter_a, bm, BM_FACES_OF_MESH) {
    BMVert *v1, *v2;
    v1 = nullptr;
    v2 = nullptr;
    BM_ITER_ELEM (v, &iter_b, f, BM_VERTS_OF_FACE) {
      if (BM_vert_edge_count(v) == 3) {
        if (v1) {
          v2 = v;
        }
        else {
          v1 = v;
        }
      }
    }
    if (v1 && v2 && (v1 != v2) && !BM_edge_exists(v1, v2)) {
      BM_face_kill(bm, f);
      BMEdge *e = BM_edge_create(bm, v1, v2, nullptr, BM_CREATE_NOP);
      BM_elem_flag_set(e, BM_ELEM_TAG, true);
    }
  }

  BM_ITER_MESH_MUTABLE (ed, ed_next, &iter_a, bm, BM_EDGES_OF_MESH) {
    if (BM_elem_flag_test(ed, BM_ELEM_TAG)) {
      float co[3];
      mid_v3_v3v3(co, ed->v1->co, ed->v2->co);
      BMVert *vc = BM_edge_collapse(bm, ed, ed->v1, true, true);
      copy_v3_v3(vc->co, co);
    }
  }

  /* Delete faces with a 3 edge pole in all their vertices */
  BM_mesh_elem_hflag_disable_all(bm, BM_VERT | BM_EDGE | BM_FACE, BM_ELEM_TAG, false);
  BM_ITER_MESH (f, &iter_a, bm, BM_FACES_OF_MESH) {
    bool dissolve = true;
    BM_ITER_ELEM (v, &iter_b, f, BM_VERTS_OF_FACE) {
      if (BM_vert_edge_count(v) != 3) {
        dissolve = false;
      }
    }
    if (dissolve) {
      BM_ITER_ELEM (v, &iter_b, f, BM_VERTS_OF_FACE) {
        BM_elem_flag_set(v, BM_ELEM_TAG, true);
      }
    }
  }
  BM_mesh_delete_hflag_context(bm, BM_ELEM_TAG, DEL_VERTS);

  BM_ITER_MESH (ed, &iter_a, bm, BM_EDGES_OF_MESH) {
    if (BM_edge_face_count(ed) != 2) {
      BM_elem_flag_set(ed, BM_ELEM_TAG, true);
    }
  }
  BM_mesh_edgenet(bm, false, true);

  /* Smooth the result */
  for (int i = 0; i < 4; i++) {
    BM_ITER_MESH (v, &iter_a, bm, BM_VERTS_OF_MESH) {
      float co[3];
      zero_v3(co);
      BM_ITER_ELEM (ed, &iter_b, v, BM_EDGES_OF_VERT) {
        BMVert *vert = BM_edge_other_vert(ed, v);
        add_v3_v3(co, vert->co);
      }
      mul_v3_fl(co, 1.0f / float(BM_vert_edge_count(v)));
      mid_v3_v3v3(v->co, v->co, co);
    }
  }

  BM_mesh_normals_update(bm);

  BM_mesh_elem_hflag_disable_all(bm, BM_VERT | BM_EDGE | BM_FACE, BM_ELEM_SELECT, false);
  BM_mesh_elem_hflag_enable_all(bm, BM_FACE, BM_ELEM_TAG, false);
  BMO_op_callf(bm,
               (BMO_FLAG_DEFAULTS & ~BMO_FLAG_RESPECT_HIDE),
               "recalc_face_normals faces=%hf",
               BM_ELEM_TAG);
  BM_mesh_elem_hflag_disable_all(bm, BM_VERT | BM_EDGE | BM_FACE, BM_ELEM_TAG, false);

  BMeshToMeshParams bmesh_to_mesh_params{};
  bmesh_to_mesh_params.calc_object_remap = false;
  Mesh *result = BKE_mesh_from_bmesh_nomain(bm, &bmesh_to_mesh_params, mesh);

  BM_mesh_free(bm);
  return result;
}