IsotropicRemeshing.h

#pragma once
 
#include <wmtk/utils/PartitionMesh.h>
#include <wmtk/utils/VectorUtils.h>
#include <wmtk/envelope/BoundaryEnvelope.hpp>
#include <wmtk/envelope/Envelope.hpp>
#include "wmtk/AttributeCollection.hpp"
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <igl/write_triangle_mesh.h>
#include <wmtk/threading/enumerable_thread_specific.hpp>
#include <fastenvelope/FastEnvelope.h>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <atomic>
#include <memory>
#include <optional>
#include <queue>
 
namespace wmtk::components::isotropic_remeshing {
 
struct VertexAttributes
{
    Eigen::Vector3d pos;
    // TODO: in fact, partition id should not be vertex attribute, it is a fixed marker to distinguish tuple/operations.
    size_t partition_id;
    bool freeze = false;
    bool input_boundary = false;
};
 
struct FeatureAttributes
{
    std::vector<size_t> nbrs;
    std::vector<int64_t> ids;
    bool fixed = false;
    int64_t vertex_id = -1;
};
 
struct FaceFieldAttributes
{
    std::vector<double> values;
    int64_t region = -1;
};
 
struct ReplacedSurface
{
    std::vector<std::pair<int64_t, Eigen::Vector3d>> normals;
    size_t against = 0;
    size_t star_folds = 0;
    struct Crease
    {
        size_t a = 0, b = 0;
        int orientation = 0;
    };
    std::vector<Crease> creases;
    std::optional<Eigen::Vector3d> normal(int64_t region) const
    {
        for (const auto& [r, n] : normals) {
            if (r == region) return n;
        }
        return std::nullopt;
    }
    const Crease* crease(size_t a, size_t b) const
    {
        for (const Crease& c : creases) {
            if ((c.a == a && c.b == b) || (c.a == b && c.b == a)) return &c;
        }
        return nullptr;
    }
    void rename(size_t from, size_t to)
    {
        for (Crease& c : creases) {
            if (c.a == from) c.a = to;
            if (c.b == from) c.b = to;
        }
    }
    bool folded() const { return against > 0 || star_folds > 0; }
};
 
class IsotropicRemeshing : public wmtk::TriMesh
{
public:
    wmtk::Envelope m_envelope;
    bool m_has_envelope = false;
    wmtk::BoundaryEnvelope m_boundary_envelope;
    wmtk::Envelope m_feature_envelope;
 
    using VertAttCol = wmtk::AttributeCollection<VertexAttributes>;
    VertAttCol vertex_attrs;
    wmtk::AttributeCollection<FeatureAttributes> feature_attrs;
    wmtk::AttributeCollection<FaceFieldAttributes> face_field_attrs;
 
    void set_feature_edges(const std::vector<std::array<size_t, 2>>& edges);
    void set_feature_edges(
        const std::vector<std::array<size_t, 2>>& edges,
        const std::vector<int64_t>& ids);
    std::vector<std::array<size_t, 2>> feature_edges() const;
    bool is_feature_edge(size_t a, size_t b) const;
    int64_t feature_edge_id(size_t a, size_t b) const;
    size_t feature_degree(size_t v) const;
 
    void set_fixed_vertices(const std::vector<std::pair<size_t, int64_t>>& vertices);
    bool is_fixed(size_t v) const { return m_has_fixed && feature_attrs.at(v).fixed; }
    bool stays_put(size_t v) const;
    bool may_collapse_away(const Tuple& t) const;
 
    void set_face_fields(const std::vector<std::vector<double>>& values);
    bool has_face_fields() const { return m_has_face_fields; }
 
    void set_feature_envelope(double eps);
    void consolidate_mesh();
    size_t remove_feature_ears();
    bool is_feature_ear_tip(const Tuple& v) const;
 
    void set_check_inversions(bool check) { m_check_inversions = check; }
    void set_swap_quality(bool quality) { m_swap_quality = quality; }
    void set_limit_collapse_length(bool limit) { m_limit_collapse_length = limit; }
    void set_cleanup_angle(double degrees) { m_cleanup_angle = degrees; }
    size_t cleanup_triangles(double degrees);
    size_t refused_inversions() const { return m_refused_inversions; }
    bool smooth_vertex_to(size_t v, const Eigen::Vector3d& p);
 
    int retry_limit = 10;
    IsotropicRemeshing(
        std::vector<Eigen::Vector3d> _m_vertex_positions,
        int num_threads = 1,
        bool use_exact = true);
 
    ~IsotropicRemeshing() {}
 
    void create_mesh(
        size_t n_vertices,
        const std::vector<std::array<size_t, 3>>& tris,
        const std::vector<size_t>& frozen_verts = std::vector<size_t>(),
        bool m_freeze = true,
        double eps = 0,
        double boundary_eps = 0);
 
    struct PositionInfoCache
    {
        Eigen::Vector3d v1p;
        Eigen::Vector3d v2p;
        size_t partition_id;
        // Collapse: of t.vid() (v1) and of the other endpoint (v2). All false without a
        // boundary envelope, where no vertex is flagged.
        bool v1_input_boundary = false;
        bool v2_input_boundary = false;
        // An input_boundary vertex that is still on the boundary.
        bool v1_on_boundary = false;
        bool v2_on_boundary = false;
        // Split: the new vertex lands on a boundary edge that touches an input_boundary vertex.
        bool split_input_boundary = false;
        // Feature curves: the edge's endpoints (t.vid() first), whether it is a feature edge and
        // its id, and whether a collapse leaves the survivor where it is.
        size_t v1 = 0, v2 = 0;
        bool feature_edge = false;
        int64_t feature_id = -1;
        bool keep_v2_position = false;
        // Collapse and swap, checking inversions: what the operation replaces.
        ReplacedSurface replaced;
        // A split is under way: its invariants skip the feature envelope.
        bool splitting = false;
        // Collapse and swap, in cleanup_triangles(): where what the operation replaces folds
        // already, its triangles that face against the input; 0 elsewhere, where it is not read.
        size_t input_folds = 0;
    };
    wmtk::threading::enumerable_thread_specific<PositionInfoCache> position_cache;
 
    void cache_edge_positions(const Tuple& t);
 
    bool invariants(const std::vector<Tuple>& new_tris) override;
 
    // TODO: this should not be here
    void partition_mesh();
 
    // TODO: morton should not be here, but inside wmtk
    void partition_mesh_morton();
 
    size_t get_partition_id(const Tuple& loc) const
    {
        return vertex_attrs[loc.vid(*this)].partition_id;
    }
 
    bool smooth_all_vertices();
 
    Eigen::Vector3d smooth(const Tuple& t);
 
 
    Eigen::Vector3d tangential_smooth(const Tuple& t);
 
    bool collapse_edge_before(const Tuple& t) override;
    bool collapse_edge_after(const Tuple& t) override;
 
    bool swap_edge_before(const Tuple& t) override;
    bool swap_edge_after(const Tuple& t) override;
 
    std::vector<TriMesh::Tuple> new_edges_after(const std::vector<TriMesh::Tuple>& tris) const;
    std::vector<TriMesh::Tuple> new_edges_after_swap(const TriMesh::Tuple& t) const;
    std::vector<TriMesh::Tuple> replace_edges_after_split(
        const std::vector<TriMesh::Tuple>& tris,
        const size_t vid_threshold) const;
    std::vector<TriMesh::Tuple> new_sub_edges_after_split(
        const std::vector<TriMesh::Tuple>& tris) const;
 
 
    bool split_edge_before(const Tuple& t) override;
    bool split_edge_after(const Tuple& t) override;
 
    bool smooth_before(const Tuple& t) override;
    bool smooth_after(const Tuple& t) override;
 
    double compute_edge_cost_collapse(const TriMesh::Tuple& t, double L) const;
    double compute_edge_cost_split(const TriMesh::Tuple& t, double L) const;
    double compute_vertex_valence(const TriMesh::Tuple& t) const;
    double compute_swap_quality(const TriMesh::Tuple& t) const;
    std::vector<double> average_len_valen();
    bool split_remeshing(double L);
    bool collapse_remeshing(double L);
    bool swap_remeshing();
    bool uniform_remeshing(double L, int interations);
    bool write_triangle_mesh(std::string path);
 
private:
    wmtk::AttributeContainerGroup m_vertex_group;
    bool m_has_features = false;
    bool m_has_fixed = false;
    bool m_has_edge_ids = false;
    bool m_has_face_fields = false;
    bool m_has_feature_envelope = false;
    bool m_removing_ears = false;
    bool m_check_inversions = true;
    bool m_swap_quality = false;
    bool m_limit_collapse_length = false;
    double m_cleanup_angle = 0;
    double m_cleanup_floor = 0;
    bool m_cleaning = false;
    std::optional<Eigen::Vector3d> m_smooth_target;
    double m_envelope_eps = 0;
    std::vector<Eigen::Vector3d> m_input_area;
    std::vector<std::vector<double>> m_input_values;
    double m_collapse_target = 0;
    std::atomic<size_t> m_refused_inversions = 0;
 
    bool m_has_regions = false;
    void set_regions();
    int64_t region(const Tuple& f) const;
    std::vector<std::pair<int64_t, std::vector<Tuple>>> by_region(
        const std::vector<Tuple>& tris) const;
    template <typename F>
    void for_each_region(const std::vector<Tuple>& tris, F&& side) const;
    std::vector<std::array<size_t, 6>> feature_wedges(const std::vector<Tuple>& tris) const;
    bool closes_a_crease(const std::vector<Tuple>& tris, const ReplacedSurface& before) const;
    Eigen::Vector3d area_normal(const std::vector<Tuple>& tris) const;
    size_t not_facing(const std::vector<Tuple>& tris, const Eigen::Vector3d& n) const;
    size_t star_folds(const std::vector<Tuple>& tris) const;
    ReplacedSurface replaced(const std::vector<Tuple>& tris) const;
    bool more_folds(const std::vector<Tuple>& tris, const ReplacedSurface& before) const;
    bool folds(const std::vector<Tuple>& tris, const ReplacedSurface& before);
 
    double min_angle(const std::vector<Tuple>& tris) const;
    Eigen::Vector3d input_normal(const Tuple& f) const;
    size_t input_folds(const std::vector<Tuple>& tris) const;
    bool cleanup_allows(
        const std::vector<Tuple>& tris,
        const ReplacedSurface& before,
        size_t input_folds_before,
        double floor) const;
    bool cleanup_move(size_t v, double floor);
};
 
} // namespace wmtk::components::isotropic_remeshing