added STL file support, deduplication algo

Readme.md

@@ -30,4 +30,5 @@ Best used with glm and glow.
 * paired_with function for smart range
 * operator +-*/ for attributes (similar to valarray)
 * dual mesh construction
-* cast<> function
\ No newline at end of file
+* cast<> function
+* surface tracing
\ No newline at end of file

src/polymesh/algorithms.hh

@@ -17,7 +17,7 @@
 // Basic mesh operations, including:
 // - elementary subdivision
 // - intersections
-#include "algorithms/operations.hh"
+#include "algorithms/remeshing/triangulate.hh"
 
 // Mesh statistics
 #include "algorithms/components.hh"

src/polymesh/algorithms/operations.hh

 #pragma once
 
-#include "../Mesh.hh"
-
-#include <glm/glm.hpp>
+#include <polymesh/Mesh.hh>
 
 // Basic mesh operations, including:
 // - elementary subdivision
@@ -10,26 +8,43 @@
 
 namespace polymesh
 {
-/// Given a flat polymesh with convex faces, naively triangulates all faces
-void triangulate_naive(Mesh& m);
+/// Removes all faces of a given mesh
+/// NOTE: does NOT compactify!
+void remove_faces(Mesh& m);
+
+/// Removes all edges and faces of a given mesh
+/// NOTE: does NOT compactify!
+void remove_edges_and_faces(Mesh& m);
 
 /// ======== IMPLEMENTATION ========
 
-inline void triangulate_naive(Mesh& m)
+inline void remove_faces(Mesh& m)
 {
-    std::vector<vertex_handle> vs;
+    auto ll = low_level_api(m);
+
+    // set faces to removed
     for (auto f : m.faces())
-    {
-        vs = f.vertices().to_vector();
-        if (vs.size() <= 3)
-            continue;
-
-        // remove
-        m.faces().remove(f);
-
-        // triangulate
-        for (auto i = 2u; i < vs.size(); ++i)
-            m.faces().add(vs[0], vs[1], vs[i]);
-    }
+        ll.set_removed(f);
+
+    // remove all faces from half-edges
+    for (auto h : m.halfedges())
+        ll.face_of(h) = face_index::invalid();
+}
+
+inline void remove_edges_and_faces(Mesh& m)
+{
+    auto ll = low_level_api(m);
+
+    // set faces to removed
+    for (auto f : m.faces())
+        ll.set_removed(f);
+
+    // set edges to removed
+    for (auto e : m.edges())
+        ll.set_removed(e);
+
+    // remove all halfedges from vertices
+    for (auto v : m.vertices())
+        ll.outgoing_halfedge_of(v) = halfedge_index::invalid();
 }
 }

src/polymesh/algorithms/remeshing/triangulate.hh

+#pragma once
+
+#include <polymesh/Mesh.hh>
+
+// Basic mesh operations, including:
+// - elementary subdivision
+// - intersections
+
+namespace polymesh
+{
+/// Given a flat polymesh with convex faces, naively triangulates all faces
+void triangulate_naive(Mesh& m);
+
+/// ======== IMPLEMENTATION ========
+
+inline void triangulate_naive(Mesh& m)
+{
+    std::vector<vertex_handle> vs;
+    for (auto f : m.faces())
+    {
+        vs = f.vertices().to_vector();
+        if (vs.size() <= 3)
+            continue;
+
+        // remove
+        m.faces().remove(f);
+
+        // triangulate
+        for (auto i = 2u; i < vs.size(); ++i)
+            m.faces().add(vs[0], vs[1], vs[i]);
+    }
+}
+}

src/polymesh/algorithms/repair/deduplicate.hh

+#pragma once
+
+#include <unordered_map>
+
+#include <polymesh/Mesh.hh>
+
+#include "../operations.hh"
+
+namespace polymesh
+{
+/// Merges vertices that report the same key
+///
+/// Example usage:
+///     Mesh m;
+///     vertex_attribute<glm::vec3> pos;
+///     load_stl(file, m, pos);
+///     deduplicate(m, pos);
+///
+/// Note:
+///     preserves (first) vertex and face attributes ONLY!
+///     edge/halfedge attributes are UNDEFINED (will probably contain uncorrelated old data)
+///
+/// returns number of removed vertices (-1 if deduplication failed (e.g. due to non-manifoldness))
+template <class KeyF>
+int deduplicate(Mesh& m, KeyF&& kf);
+
+/// ======== IMPLEMENTATION ========
+
+template <class KeyF>
+int deduplicate(Mesh& m, KeyF&& kf)
+{
+    using KeyT = typename std::decay<decltype(kf(m.vertices().first()))>::type;
+
+    std::unordered_map<KeyT, vertex_index> remap;
+    auto new_idx = m.vertices().make_attribute<vertex_index>();
+
+    // calculate remapped vertices
+    for (auto v : m.vertices())
+    {
+        auto const& k = kf(v);
+        if (!remap.count(k))
+            remap[k] = v;
+
+        new_idx[v] = remap[k];
+    }
+
+    // calc face remapping
+    std::vector<vertex_index> poly_verts;
+    struct poly
+    {
+        face_index f;
+        int start;
+        int count;
+    };
+    std::vector<poly> polys;
+    polys.reserve(m.faces().size());
+    for (auto f : m.faces())
+    {
+        auto s = (int)poly_verts.size();
+        for (auto v : f.vertices())
+            poly_verts.push_back(new_idx[v]);
+        auto e = (int)poly_verts.size();
+        polys.push_back({f, s, e - s});
+    }
+
+    auto ll = low_level_api(m);
+
+    // remove everything except vertices
+    remove_edges_and_faces(m);
+
+    // clear edge vector (new edges are allocated from idx 0)
+    ll.clear_removed_edge_vector();
+
+    // add remapped faces
+    auto manifold = true;
+    for (auto const& p : polys)
+    {
+        if (ll.can_add_face(poly_verts.data() + p.start, p.count))
+            ll.add_face(poly_verts.data() + p.start, p.count, p.f);
+        else
+            manifold = false;
+    }
+
+    // remove duplicated vertices
+    int removed = 0;
+    for (auto v : m.vertices())
+        if (new_idx[v] != v)
+        {
+            m.vertices().remove(v);
+            ++removed;
+        }
+
+    return manifold ? removed : -1;
+}
+}

src/polymesh/formats/stl.cc

+#include "stl.hh"
+
+#include <cstddef>
+
+#include <fstream>
+#include <sstream>
+
+/*
+    UINT8[80] – Header
+    UINT32 – Number of triangles
+
+    foreach triangle
+        REAL32[3] – Normal vector
+        REAL32[3] – Vertex 1
+        REAL32[3] – Vertex 2
+        REAL32[3] – Vertex 3
+        UINT16 – Attribute byte count
+    end
+ */
+
+namespace polymesh
+{
+void write_stl_binary(const std::string &filename, const Mesh &mesh, const vertex_attribute<glm::vec3> &position, face_attribute<glm::vec3> const *normals)
+{
+    std::ofstream file(filename, std::ios_base::binary);
+    write_stl_binary(file, mesh, position, normals);
+}
+void write_stl_binary(std::ostream &out, const Mesh &mesh, const vertex_attribute<glm::vec3> &position, face_attribute<glm::vec3> const *normals)
+{
+    char header[80] = {};
+    uint32_t n_triangles = mesh.faces().size();
+
+    out.write(header, sizeof(header));
+    out.write((char const *)&n_triangles, sizeof(n_triangles));
+
+    for (auto f : mesh.faces())
+    {
+        auto n = f[normals];
+        out.write((char const *)&n, sizeof(n));
+
+        auto cnt = 0;
+        for (auto v : f.vertices())
+        {
+            if (cnt >= 3)
+            {
+                std::cerr << "STL only supports triangles" << std::endl;
+                break;
+            }
+
+            auto p = position[v];
+            out.write((char const *)&p, sizeof(p));
+
+            ++cnt;
+        }
+
+        uint16_t attr_cnt = 0;
+        out.write((char const *)&attr_cnt, sizeof(attr_cnt));
+    }
+}
+
+bool read_stl(const std::string &filename, Mesh &mesh, vertex_attribute<glm::vec3> &position, face_attribute<glm::vec3> *normals)
+{
+    std::ifstream file(filename);
+    if (!file.good())
+        return false;
+
+    return read_stl(file, mesh, position, normals);
+}
+
+bool read_stl(std::istream &input, Mesh &mesh, vertex_attribute<glm::vec3> &position, face_attribute<glm::vec3> *normals)
+{
+    return is_ascii_stl(input) ? read_stl_ascii(input, mesh, position, normals) : read_stl_binary(input, mesh, position, normals);
+}
+
+bool read_stl_binary(std::istream &input, Mesh &mesh, vertex_attribute<glm::vec3> &position, face_attribute<glm::vec3> *normals)
+{
+    mesh.clear();
+
+    char header[80];
+    input.read(header, sizeof(header));
+
+    // if (header[0] == 's' && header[1] == 'o' && header[2] == 'l' && header[3] == 'i' && header[4] == 'd')
+    // {
+    //     std::cerr << "ASCII STL is not supported" << std::endl;
+    //     return false;
+    // }
+
+    uint32_t n_triangles;
+    input.read((char *)&n_triangles, sizeof(n_triangles));
+
+    mesh.faces().reserve(n_triangles);
+    mesh.vertices().reserve(n_triangles * 3);
+    mesh.halfedges().reserve(n_triangles * 3);
+
+    for (auto i = 0; i < n_triangles; ++i)
+    {
+        if (!input.good())
+        {
+            std::cerr << "Premature end of file";
+            return false;
+        }
+
+        auto v0 = mesh.vertices().add();
+        auto v1 = mesh.vertices().add();
+        auto v2 = mesh.vertices().add();
+        auto f = mesh.faces().add(v0, v1, v2);
+
+        input.read((char *)&f[normals], sizeof(glm::vec3));
+        input.read((char *)&position[v0], sizeof(glm::vec3));
+        input.read((char *)&position[v1], sizeof(glm::vec3));
+        input.read((char *)&position[v2], sizeof(glm::vec3));
+        uint16_t attr_cnt;
+        input.read((char *)&attr_cnt, sizeof(attr_cnt));
+    }
+
+    return true;
+}
+
+bool read_stl_ascii(std::istream &input, Mesh &mesh, vertex_attribute<glm::vec3> &position, face_attribute<glm::vec3> *normals)
+{
+    mesh.clear();
+
+    std::string s;
+    input >> s;
+
+    if (s != "solid")
+    {
+        std::cerr << "File does not seem to be ASCII stl" << std::endl;
+        return false;
+    }
+
+    // name
+    do
+    {
+        input >> s;
+    } while (input.good() && s != "endsolid" && s != "facet");
+
+    while (input.good() && s != "endsolid")
+    {
+        assert(s == "facet");
+
+        vertex_handle v[3];
+        v[0] = mesh.vertices().add();
+        v[1] = mesh.vertices().add();
+        v[2] = mesh.vertices().add();
+        auto f = mesh.faces().add(v);
+
+        input >> s;
+        assert(s == "normal");
+        glm::vec3 n;
+        input >> n.x >> n.y >> n.z;
+        f[normals] = n;
+
+        input >> s;
+        assert(s == "outer");
+        input >> s;
+        assert(s == "loop");
+
+        for (auto i = 0; i < 3; ++i)
+        {
+            input >> s;
+            assert(s == "vertex");
+            glm::vec3 p;
+            input >> p.x >> p.y >> p.z;
+            position[v[i]] = p;
+        }
+
+        input >> s;
+        assert(s == "endloop");
+
+        input >> s;
+        assert(s == "endfacet");
+
+        input >> s; // for next iteration
+    }
+
+    return true;
+}
+
+bool is_ascii_stl(std::istream &input)
+{
+    auto savp = input.tellg();
+
+    char solid[6];
+    input.read(solid, sizeof(solid));
+
+    if (solid[0] != 's' || solid[1] != 'o' || solid[2] != 'l' || solid[3] != 'i' || solid[4] != 'd' || solid[5] != ' ')
+    {
+        input.seekg(savp, std::ios_base::beg);
+        return false;
+    }
+
+    std::string s;
+    input >> s;
+
+    if (s == "facet" || s == "endsolid")
+    {
+        input.seekg(savp, std::ios_base::beg);
+        return true;
+    }
+
+    input >> s;
+
+    if (s == "facet" || s == "endsolid")
+    {
+        input.seekg(savp, std::ios_base::beg);
+        return true;
+    }
+
+    input.seekg(savp, std::ios_base::beg);
+    return false;
+}
+}

src/polymesh/formats/stl.hh

+#pragma once
+
+#include <glm/glm.hpp>
+
+#include <iostream>
+#include <string>
+
+#include "../Mesh.hh"
+
+namespace polymesh
+{
+void write_stl_binary(std::string const& filename, Mesh const& mesh, vertex_attribute<glm::vec3> const& position, face_attribute<glm::vec3> const* normals = nullptr);
+void write_stl_binary(std::ostream& out, Mesh const& mesh, vertex_attribute<glm::vec3> const& position, face_attribute<glm::vec3> const* normals = nullptr);
+bool read_stl(std::string const& filename, Mesh& mesh, vertex_attribute<glm::vec3>& position, face_attribute<glm::vec3>* normals = nullptr);
+bool read_stl(std::istream& input, Mesh& mesh, vertex_attribute<glm::vec3>& position, face_attribute<glm::vec3>* normals = nullptr);
+bool read_stl_binary(std::istream& input, Mesh& mesh, vertex_attribute<glm::vec3>& position, face_attribute<glm::vec3>* normals = nullptr);
+bool read_stl_ascii(std::istream& input, Mesh& mesh, vertex_attribute<glm::vec3>& position, face_attribute<glm::vec3>* normals = nullptr);
+bool is_ascii_stl(std::istream& input);
+}

src/polymesh/impl/impl_low_level_api_base.hh

@@ -89,6 +89,30 @@ int low_level_api_base<MeshT>::size_valid_halfedges() const
     return m.size_valid_halfedges();
 }
 
+template <class MeshT>
+int low_level_api_base<MeshT>::size_removed_faces() const
+{
+    return m.size_all_faces() - m.size_removed_faces();
+}
+
+template <class MeshT>
+int low_level_api_base<MeshT>::size_removed_vertices() const
+{
+    return m.size_all_vertices() - m.size_removed_vertices();
+}
+
+template <class MeshT>
+int low_level_api_base<MeshT>::size_removed_edges() const
+{
+    return m.size_all_edges() - m.size_removed_edges();
+}
+
+template <class MeshT>
+int low_level_api_base<MeshT>::size_removed_halfedges() const
+{
+    return m.size_all_halfedges() - m.size_removed_halfedges();
+}
+
 template <class MeshT>
 bool low_level_api_base<MeshT>::can_add_face(const vertex_handle *v_handles, int vcnt) const
 {

src/polymesh/impl/impl_low_level_api_mutable.hh

@@ -20,35 +20,43 @@ inline void low_level_api_mutable::permute_faces(const std::vector<int> &p) cons
 inline void low_level_api_mutable::permute_edges(const std::vector<int> &p) const { m.permute_edges(p); }
 inline void low_level_api_mutable::permute_vertices(const std::vector<int> &p) const { m.permute_vertices(p); }
 
-inline face_index low_level_api_mutable::add_face(const vertex_handle *v_handles, int vcnt) const
+inline face_index low_level_api_mutable::add_face(const vertex_handle *v_handles, int vcnt, face_index res_idx) const
 {
     m.mFaceInsertCache.resize(vcnt);
     for (auto i = 0; i < vcnt; ++i)
         m.mFaceInsertCache[i] = add_or_get_halfedge(v_handles[i].idx, v_handles[(i + 1) % vcnt].idx);
-    return add_face(m.mFaceInsertCache.data(), vcnt);
+    return add_face(m.mFaceInsertCache.data(), vcnt, res_idx);
 }
 
-inline face_index low_level_api_mutable::add_face(const vertex_index *v_indices, int vcnt) const
+inline face_index low_level_api_mutable::add_face(const vertex_index *v_indices, int vcnt, face_index res_idx) const
 {
     m.mFaceInsertCache.resize(vcnt);
     for (auto i = 0; i < vcnt; ++i)
         m.mFaceInsertCache[i] = add_or_get_halfedge(v_indices[i], v_indices[(i + 1) % vcnt]);
-    return add_face(m.mFaceInsertCache.data(), vcnt);
+    return add_face(m.mFaceInsertCache.data(), vcnt, res_idx);
 }
 
-inline face_index low_level_api_mutable::add_face(const halfedge_handle *half_loop, int vcnt) const
+inline face_index low_level_api_mutable::add_face(const halfedge_handle *half_loop, int vcnt, face_index res_idx) const
 {
     m.mFaceInsertCache.resize(vcnt);
     for (auto i = 0; i < vcnt; ++i)
         m.mFaceInsertCache[i] = half_loop[i].idx;
-    return add_face(m.mFaceInsertCache.data(), vcnt);
+    return add_face(m.mFaceInsertCache.data(), vcnt, res_idx);
 }
 
-inline face_index low_level_api_mutable::add_face(const halfedge_index *half_loop, int vcnt) const
+inline face_index low_level_api_mutable::add_face(const halfedge_index *half_loop, int vcnt, face_index res_idx) const
 {
     assert(vcnt >= 3 && "no support for less-than-triangular faces");
+    assert(res_idx.is_invalid() || is_removed(res_idx) && "resurrected index must be previously removed!");
 
-    auto fidx = alloc_face();
+    auto fidx = res_idx.is_valid() ? res_idx : alloc_face();
+
+    // on resurrect: fix counts
+    if (res_idx.is_valid())
+    {
+        m.mRemovedFaces--;
+        // no mCompact change!
+    }
 
     // ensure that half-edges are adjacent at each vertex
     for (auto i = 0; i < vcnt; ++i)
@@ -322,6 +330,19 @@ inline void low_level_api_mutable::remove_vertex(vertex_index v_idx) const
     set_removed(v_idx);
 }
 
+inline void low_level_api_mutable::clear_removed_edge_vector() const
+{
+    assert(m.edges().empty() && "only works for no-edge meshes");
+
+    m.mHalfedgeToFace.clear();
+    m.mHalfedgeToNextHalfedge.clear();
+    m.mHalfedgeToPrevHalfedge.clear();
+    m.mHalfedgeToVertex.clear();
+
+    m.mRemovedHalfedges = 0;
+    // no mCompact change!
+}
+
 inline void low_level_api_mutable::fix_boundary_state_of(vertex_index v_idx) const
 {
     assert(!is_isolated(v_idx));

src/polymesh/low_level_api.hh

@@ -50,6 +50,11 @@ public:
     int size_valid_edges() const;
     int size_valid_halfedges() const;
 
+    int size_removed_faces() const;
+    int size_removed_vertices() const;
+    int size_removed_edges() const;
+    int size_removed_halfedges() const;
+
     // traversal helper
 public:
     // returns the next valid idx (returns the given one if valid)
@@ -164,10 +169,11 @@ public:
     /// Adds a face consisting of N vertices
     /// The vertices must already be sorted in CCW order
     /// (note: trying to add already existing halfedges triggers assertions)
-    face_index add_face(vertex_handle const* v_handles, int vcnt) const;
-    face_index add_face(vertex_index const* v_indices, int vcnt) const;
-    face_index add_face(halfedge_handle const* half_loop, int vcnt) const;
-    face_index add_face(halfedge_index const* half_loop, int vcnt) const;
+    /// Optional: provide an index that should be "resurrected" (must be a removed slot)
+    face_index add_face(vertex_handle const* v_handles, int vcnt, face_index res_idx = {}) const;
+    face_index add_face(vertex_index const* v_indices, int vcnt, face_index res_idx = {}) const;
+    face_index add_face(halfedge_handle const* half_loop, int vcnt, face_index res_idx = {}) const;
+    face_index add_face(halfedge_index const* half_loop, int vcnt, face_index res_idx = {}) const;
 
     /// Adds an edge between two existing, distinct vertices
     /// if edge already exists, returns it
@@ -207,6 +213,11 @@ public:
     /// removes all adjacent edges, then the vertex
     void remove_vertex(vertex_index v_idx) const;
 
+    /// special purpose function:
+    /// CAUTION: only works if edges.size() == 0
+    /// clears the edge vector
+    void clear_removed_edge_vector() const;
+
     /// Overrides the saved number of removed primitives
     /// CAUTION: only use if you know what you do!
     void set_removed_counts(int r_vertices, int r_faces, int r_edges);