test/source/blender/blenkernel/intern/mesh_fair.cc

/* SPDX-License-Identifier: GPL-2.0-or-later */

/** \file
 * \ingroup bke
 * Mesh Fairing algorithm designed by Brett Fedack, used in the addon "Mesh Fairing":
 * https://github.com/fedackb/mesh-fairing.
 */

#include "BLI_map.hh"
#include "BLI_math.h"
#include "BLI_vector.hh"

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

#include "BKE_lib_id.h"
#include "BKE_lib_query.h"
#include "BKE_mesh.h"
#include "BKE_mesh_fair.h"
#include "BKE_mesh_mapping.h"

#include "bmesh.h"
#include "bmesh_tools.h"

#include "MEM_guardedalloc.h"
#include "eigen_capi.h"

using blender::Map;
using blender::MutableSpan;
using blender::Span;
using blender::Vector;
using std::array;

class VertexWeight {
 public:
  virtual float weight_at_index(const int index) = 0;
  virtual ~VertexWeight() = default;
};

class LoopWeight {
 public:
  virtual float weight_at_index(const int index) = 0;
  virtual ~LoopWeight() = default;
};

class FairingContext {
 public:
  /* Get coordinates of vertices which are adjacent to the loop with specified index. */
  virtual void adjacents_coords_from_loop(const int loop,
                                          float r_adj_next[3],
                                          float r_adj_prev[3]) = 0;

  /* Get the other vertex index for a loop. */
  virtual int other_vertex_index_from_loop(const int loop, const uint v) = 0;

  int vertex_count_get()
  {
    return totvert_;
  }

  int loop_count_get()
  {
    return totvert_;
  }

  MeshElemMap *vertex_loop_map_get(const int v)
  {
    return &vlmap_[v];
  }

  float *vertex_deformation_co_get(const int v)
  {
    return co_[v];
  }

  virtual ~FairingContext() = default;

  void fair_vertices(bool *affected,
                     const eMeshFairingDepth depth,
                     VertexWeight *vertex_weight,
                     LoopWeight *loop_weight)
  {

    fair_vertices_ex(affected, (int)depth, vertex_weight, loop_weight);
  }

 protected:
  Vector<float *> co_;

  int totvert_;
  int totloop_;

  MeshElemMap *vlmap_;
  int *vlmap_mem_;

 private:
  void fair_setup_fairing(const int v,
                          const int i,
                          LinearSolver *solver,
                          float multiplier,
                          const int depth,
                          Map<int, int> &vert_col_map,
                          VertexWeight *vertex_weight,
                          LoopWeight *loop_weight)
  {
    if (depth == 0) {
      if (vert_col_map.contains(v)) {
        const int j = vert_col_map.lookup(v);
        EIG_linear_solver_matrix_add(solver, i, j, -multiplier);
        return;
      }
      for (int j = 0; j < 3; j++) {
        EIG_linear_solver_right_hand_side_add(solver, j, i, multiplier * co_[v][j]);
      }
      return;
    }

    float w_ij_sum = 0;
    const float w_i = vertex_weight->weight_at_index(v);
    MeshElemMap *vlmap_elem = &vlmap_[v];
    for (int l = 0; l < vlmap_elem->count; l++) {
      const int l_index = vlmap_elem->indices[l];
      const int other_vert = other_vertex_index_from_loop(l_index, v);
      const float w_ij = loop_weight->weight_at_index(l_index);
      w_ij_sum += w_ij;
      fair_setup_fairing(other_vert,
                         i,
                         solver,
                         w_i * w_ij * multiplier,
                         depth - 1,
                         vert_col_map,
                         vertex_weight,
                         loop_weight);
    }
    fair_setup_fairing(v,
                       i,
                       solver,
                       -1 * w_i * w_ij_sum * multiplier,
                       depth - 1,
                       vert_col_map,
                       vertex_weight,
                       loop_weight);
  }

  void fair_vertices_ex(const bool *affected,
                        const int order,
                        VertexWeight *vertex_weight,
                        LoopWeight *loop_weight)
  {
    Map<int, int> vert_col_map;
    int num_affected_vertices = 0;
    for (int i = 0; i < totvert_; i++) {
      if (!affected[i]) {
        continue;
      }
      vert_col_map.add(i, num_affected_vertices);
      num_affected_vertices++;
    }

    /* Early return, nothing to do. */
    if (ELEM(num_affected_vertices, 0, totvert_)) {
      return;
    }

    /* Setup fairing matrices */
    LinearSolver *solver = EIG_linear_solver_new(num_affected_vertices, num_affected_vertices, 3);
    for (auto item : vert_col_map.items()) {
      const int v = item.key;
      const int col = item.value;
      fair_setup_fairing(v, col, solver, 1.0f, order, vert_col_map, vertex_weight, loop_weight);
    }

    /* Solve linear system */
    EIG_linear_solver_solve(solver);

    /* Copy the result back to the mesh */
    for (auto item : vert_col_map.items()) {
      const int v = item.key;
      const int col = item.value;
      for (int j = 0; j < 3; j++) {
        co_[v][j] = EIG_linear_solver_variable_get(solver, j, col);
      }
    }

    /* Free solver data */
    EIG_linear_solver_delete(solver);
  }
};

class MeshFairingContext : public FairingContext {
 public:
  MeshFairingContext(Mesh *mesh, MVert *deform_mverts)
  {
    totvert_ = mesh->totvert;
    totloop_ = mesh->totloop;

    MutableSpan<MVert> verts = mesh->vertices_for_write();
    medge_ = mesh->edges();
    mpoly_ = mesh->polygons();
    mloop_ = mesh->loops();
    BKE_mesh_vert_loop_map_create(&vlmap_,
                                  &vlmap_mem_,
                                  mpoly_.data(),
                                  mloop_.data(),
                                  mesh->totvert,
                                  mesh->totpoly,
                                  mesh->totloop);

    /* Deformation coords. */
    co_.reserve(mesh->totvert);
    if (deform_mverts) {
      for (int i = 0; i < mesh->totvert; i++) {
        co_[i] = deform_mverts[i].co;
      }
    }
    else {
      for (int i = 0; i < mesh->totvert; i++) {
        co_[i] = verts[i].co;
      }
    }

    loop_to_poly_map_.reserve(mesh->totloop);
    for (int i = 0; i < mesh->totpoly; i++) {
      for (int l = 0; l < mpoly_[i].totloop; l++) {
        loop_to_poly_map_[l + mpoly_[i].loopstart] = i;
      }
    }
  }

  ~MeshFairingContext() override
  {
    MEM_SAFE_FREE(vlmap_);
    MEM_SAFE_FREE(vlmap_mem_);
  }

  void adjacents_coords_from_loop(const int loop,
                                  float r_adj_next[3],
                                  float r_adj_prev[3]) override
  {
    const int vert = mloop_[loop].v;
    const MPoly *p = &mpoly_[loop_to_poly_map_[loop]];
    const int corner = poly_find_loop_from_vert(p, &mloop_[p->loopstart], vert);
    copy_v3_v3(r_adj_next, co_[ME_POLY_LOOP_NEXT(mloop_, p, corner)->v]);
    copy_v3_v3(r_adj_prev, co_[ME_POLY_LOOP_PREV(mloop_, p, corner)->v]);
  }

  int other_vertex_index_from_loop(const int loop, const uint v) override
  {
    const MEdge *e = &medge_[mloop_[loop].e];
    if (e->v1 == v) {
      return e->v2;
    }
    return e->v1;
  }

 protected:
  Mesh *mesh_;
  Span<MLoop> mloop_;
  Span<MPoly> mpoly_;
  Span<MEdge> medge_;
  Vector<int> loop_to_poly_map_;
};

class BMeshFairingContext : public FairingContext {
 public:
  BMeshFairingContext(BMesh *bm)
  {
    this->bm = bm;
    totvert_ = bm->totvert;
    totloop_ = bm->totloop;

    BM_mesh_elem_table_ensure(bm, BM_VERT);
    BM_mesh_elem_index_ensure(bm, BM_LOOP);

    /* Deformation coords. */
    co_.reserve(bm->totvert);
    for (int i = 0; i < bm->totvert; i++) {
      BMVert *v = BM_vert_at_index(bm, i);
      co_[i] = v->co;
    }

    bmloop_.reserve(bm->totloop);
    vlmap_ = (MeshElemMap *)MEM_calloc_arrayN(bm->totvert, sizeof(MeshElemMap), "bmesh loop map");
    vlmap_mem_ = (int *)MEM_malloc_arrayN(bm->totloop, sizeof(int), "bmesh loop map mempool");

    BMVert *v;
    BMLoop *l;
    BMIter iter;
    BMIter loop_iter;
    int index_iter = 0;

    /* This initializes both the bmloop and the vlmap for bmesh in a single loop. */
    BM_ITER_MESH (v, &iter, bm, BM_VERTS_OF_MESH) {
      int loop_count = 0;
      const int vert_index = BM_elem_index_get(v);
      vlmap_[vert_index].indices = &vlmap_mem_[index_iter];
      BM_ITER_ELEM (l, &loop_iter, v, BM_LOOPS_OF_VERT) {
        const int loop_index = BM_elem_index_get(l);
        bmloop_[loop_index] = l;
        vlmap_mem_[index_iter] = loop_index;
        index_iter++;
        loop_count++;
      }
      vlmap_[vert_index].count = loop_count;
    }
  }

  ~BMeshFairingContext() override
  {
    MEM_SAFE_FREE(vlmap_);
    MEM_SAFE_FREE(vlmap_mem_);
  }

  void adjacents_coords_from_loop(const int loop,
                                  float r_adj_next[3],
                                  float r_adj_prev[3]) override
  {
    copy_v3_v3(r_adj_next, bmloop_[loop]->next->v->co);
    copy_v3_v3(r_adj_prev, bmloop_[loop]->prev->v->co);
  }

  int other_vertex_index_from_loop(const int loop, const uint v) override
  {
    BMLoop *l = bmloop_[loop];
    BMVert *bmvert = BM_vert_at_index(bm, v);
    BMVert *bm_other_vert = BM_edge_other_vert(l->e, bmvert);
    return BM_elem_index_get(bm_other_vert);
  }

 protected:
  BMesh *bm;
  Vector<BMLoop *> bmloop_;
};

class UniformVertexWeight : public VertexWeight {
 public:
  UniformVertexWeight(FairingContext *fairing_context)
  {
    const int totvert = fairing_context->vertex_count_get();
    vertex_weights_.reserve(totvert);
    for (int i = 0; i < totvert; i++) {
      const int tot_loop = fairing_context->vertex_loop_map_get(i)->count;
      if (tot_loop != 0) {
        vertex_weights_[i] = 1.0f / tot_loop;
      }
      else {
        vertex_weights_[i] = FLT_MAX;
      }
    }
  }

  float weight_at_index(const int index) override
  {
    return vertex_weights_[index];
  }

 private:
  Vector<float> vertex_weights_;
};

class VoronoiVertexWeight : public VertexWeight {

 public:
  VoronoiVertexWeight(FairingContext *fairing_context)
  {

    const int totvert = fairing_context->vertex_count_get();
    vertex_weights_.reserve(totvert);
    for (int i = 0; i < totvert; i++) {

      float area = 0.0f;
      float a[3];
      copy_v3_v3(a, fairing_context->vertex_deformation_co_get(i));
      const float acute_threshold = M_PI_2;

      MeshElemMap *vlmap_elem = fairing_context->vertex_loop_map_get(i);
      for (int l = 0; l < vlmap_elem->count; l++) {
        const int l_index = vlmap_elem->indices[l];

        float b[3], c[3], d[3];
        fairing_context->adjacents_coords_from_loop(l_index, b, c);

        if (angle_v3v3v3(c, fairing_context->vertex_deformation_co_get(i), b) < acute_threshold) {
          calc_circumcenter(d, a, b, c);
        }
        else {
          add_v3_v3v3(d, b, c);
          mul_v3_fl(d, 0.5f);
        }

        float t[3];
        add_v3_v3v3(t, a, b);
        mul_v3_fl(t, 0.5f);
        area += area_tri_v3(a, t, d);

        add_v3_v3v3(t, a, c);
        mul_v3_fl(t, 0.5f);
        area += area_tri_v3(a, d, t);
      }

      vertex_weights_[i] = area != 0.0f ? 1.0f / area : 1e12;
    }
  }

  float weight_at_index(const int index) override
  {
    return vertex_weights_[index];
  }

 private:
  Vector<float> vertex_weights_;

  void calc_circumcenter(float r[3], const float a[3], const float b[3], const float c[3])
  {
    float ab[3];
    sub_v3_v3v3(ab, b, a);

    float ac[3];
    sub_v3_v3v3(ac, c, a);

    float ab_cross_ac[3];
    cross_v3_v3v3(ab_cross_ac, ab, ac);

    if (len_squared_v3(ab_cross_ac) > 0.0f) {
      float d[3];
      cross_v3_v3v3(d, ab_cross_ac, ab);
      mul_v3_fl(d, len_squared_v3(ac));

      float t[3];
      cross_v3_v3v3(t, ac, ab_cross_ac);
      mul_v3_fl(t, len_squared_v3(ab));

      add_v3_v3(d, t);

      mul_v3_fl(d, 1.0f / (2.0f * len_squared_v3(ab_cross_ac)));

      add_v3_v3v3(r, a, d);
      return;
    }
    copy_v3_v3(r, a);
  }
};

class UniformLoopWeight : public LoopWeight {
 public:
  float weight_at_index(const int UNUSED(index)) override
  {
    return 1.0f;
  }
};

static void prefair_and_fair_vertices(FairingContext *fairing_context,
                                      bool *affected_vertices,
                                      const eMeshFairingDepth depth)
{
  /* Prefair. */
  UniformVertexWeight *uniform_vertex_weights = new UniformVertexWeight(fairing_context);
  UniformLoopWeight *uniform_loop_weights = new UniformLoopWeight();
  fairing_context->fair_vertices(
      affected_vertices, depth, uniform_vertex_weights, uniform_loop_weights);
  delete uniform_vertex_weights;

  /* Fair. */
  VoronoiVertexWeight *voronoi_vertex_weights = new VoronoiVertexWeight(fairing_context);
  /* TODO: Implement cotangent loop weights. */
  fairing_context->fair_vertices(
      affected_vertices, depth, voronoi_vertex_weights, uniform_loop_weights);

  delete uniform_loop_weights;
  delete voronoi_vertex_weights;
}

void BKE_mesh_prefair_and_fair_vertices(struct Mesh *mesh,
                                        struct MVert *deform_mverts,
                                        bool *affect_vertices,
                                        const eMeshFairingDepth depth)
{
  MeshFairingContext *fairing_context = new MeshFairingContext(mesh, deform_mverts);
  prefair_and_fair_vertices(fairing_context, affect_vertices, depth);
  delete fairing_context;
}

void BKE_bmesh_prefair_and_fair_vertices(struct BMesh *bm,
                                         bool *affect_vertices,
                                         const eMeshFairingDepth depth)
{
  BMeshFairingContext *fairing_context = new BMeshFairingContext(bm);
  prefair_and_fair_vertices(fairing_context, affect_vertices, depth);
  delete fairing_context;
}