working on permutations and layout optimization

Readme.md

@@ -25,4 +25,6 @@ Best used with glm and glow.
 * vertex split?
 * half-edge collapse
 * normal, tangent, bitangent computation
-* attribute iterator
\ No newline at end of file
+* attribute iterator
+* primitive sort functions, better remap function, cache optimization
+* structure of arrays instead of AOS
\ No newline at end of file

src/polymesh/Mesh.hh

@@ -244,6 +244,13 @@ private:
     /// returns the vertex that this half-edge is leaving from
     vertex_index from_vertex_of(halfedge_index idx) const { return halfedge(opposite(idx)).to_vertex; }
 
+    /// applies an index remapping to all face indices (p[curr_idx] = new_idx)
+    void permute_faces(std::vector<int> const& p);
+    /// applies an index remapping to all edge (and half-edge) indices (p[curr_idx] = new_idx)
+    void permute_edges(std::vector<int> const& p);
+    /// applies an index remapping to all vertices indices (p[curr_idx] = new_idx)
+    void permute_vertices(std::vector<int> const& p);
+
     // internal datastructures
 private:
     // 4 byte per face

src/polymesh/algorithms/optimization.hh

+#pragma once
+
+#include "../Mesh.hh"
+
+#include "../detail/permutation.hh"
+#include "../detail/random.hh"
+
+namespace polymesh
+{
+/// Optimizes mesh layout for face traversals
+// TODO: half-edge iteration should be cache local
+void optimize_for_face_traversal(Mesh& m);
+
+/// Optimizes mesh layout for vertex traversals
+// TODO: half-edge iteration should be cache local
+void optimize_for_vertex_traversal(Mesh& m);
+
+/// Optimizes mesh layout for indexed face rendering
+void optimize_for_rendering(Mesh& m);
+
+/// optimizes edge indices for a given face neighborhood
+void optimize_edges_for_faces(Mesh& m);
+/// optimizes vertex indices for a given face neighborhood
+void optimize_vertices_for_faces(Mesh& m);
+
+/// optimizes edge indices for a given vertex neighborhood
+void optimize_edges_for_vertices(Mesh& m);
+/// optimizes face indices for a given vertex neighborhood
+void optimize_faces_for_vertices(Mesh& m);
+
+/// Calculates a cache-coherent face layout in O(n log n) time
+/// Can be applied using m.faces().permute(...)
+/// Returns remapping [curr_idx] = new_idx
+std::vector<int> cache_coherent_face_layout(Mesh const& m);
+
+/// Calculates a cache-coherent vertex layout in O(n log n) time
+/// Can be applied using m.vertices().permute(...)
+/// Returns remapping [curr_idx] = new_idx
+std::vector<int> cache_coherent_vertex_layout(Mesh const& m);
+
+/// ======== IMPLEMENTATION ========
+
+inline void optimize_for_face_traversal(Mesh& m)
+{
+    m.faces().permute(cache_coherent_face_layout(m));
+    optimize_edges_for_faces(m);
+    optimize_vertices_for_faces(m);
+}
+
+inline void optimize_for_vertex_traversal(Mesh& m)
+{
+    m.vertices().permute(cache_coherent_vertex_layout(m));
+    optimize_edges_for_vertices(m);
+    optimize_faces_for_vertices(m);
+}
+
+inline void optimize_for_rendering(Mesh& m)
+{
+    // TODO
+    assert(0 && "TODO");
+}
+
+inline std::vector<int> cache_coherent_face_layout(Mesh const& m)
+{
+    std::vector<int> id;
+    for (auto i = 0u; i < m.faces().size(); ++i)
+        id.push_back(i);
+    // TODO
+    return id;
+}
+
+inline std::vector<int> cache_coherent_vertex_layout(Mesh const& m)
+{
+    assert(0 && "TODO");
+    return {};
+}
+
+inline void optimize_edges_for_faces(Mesh& m)
+{
+    uint64_t rng = 1;
+
+    std::vector<std::pair<int, int>> face_edge_indices;
+    for (auto e : m.edges())
+    {
+        auto fA = e.faceA();
+        auto fB = e.faceB();
+        auto f = fA.is_invalid() ? fB : fB.is_invalid() ? fA : detail::xorshift64star(rng) % 2 ? fA : fB;
+        face_edge_indices.emplace_back(f.idx.value, e.idx.value);
+    }
+
+    // sort by face idx
+    sort(face_edge_indices.begin(), face_edge_indices.end());
+
+    // extract edge indices
+    std::vector<int> permutation(face_edge_indices.size());
+    for (auto i = 0u; i < face_edge_indices.size(); ++i)
+        permutation[i] = face_edge_indices[i].second;
+
+    // apply permutation
+    m.edges().permute(permutation);
+}
+
+inline void optimize_edges_for_vertices(Mesh& m)
+{
+    uint64_t rng = 1;
+
+    std::vector<std::pair<int, int>> vertex_edge_indices;
+    for (auto e : m.edges())
+    {
+        auto r = detail::xorshift64star(rng);
+        vertex_edge_indices.emplace_back(r % 2 ? e.vertexA().idx.value : e.vertexB().idx.value, e.idx.value);
+    }
+
+    // sort by vertex idx
+    sort(vertex_edge_indices.begin(), vertex_edge_indices.end());
+
+    // extract edge indices
+    std::vector<int> permutation(vertex_edge_indices.size());
+    for (auto i = 0u; i < vertex_edge_indices.size(); ++i)
+        permutation[i] = vertex_edge_indices[i].second;
+
+    // apply permutation
+    m.edges().permute(permutation);
+}
+
+inline void optimize_faces_for_vertices(Mesh& m)
+{
+    uint64_t rng = 1;
+
+    std::vector<std::pair<int, int>> vertex_face_indices;
+    for (auto f : m.faces())
+    {
+        vertex_handle vv;
+        auto cnt = 0;
+        for (auto v : f.vertices())
+        {
+            ++cnt;
+            if (detail::xorshift64star(rng) % cnt == 0)
+                vv = v;
+        }
+        vertex_face_indices.emplace_back(vv.idx.value, f.idx.value);
+    }
+
+    // sort by vertex idx
+    sort(vertex_face_indices.begin(), vertex_face_indices.end());
+
+    // extract face indices
+    std::vector<int> permutation(vertex_face_indices.size());
+    for (auto i = 0u; i < vertex_face_indices.size(); ++i)
+        permutation[i] = vertex_face_indices[i].second;
+
+    // apply permutation
+    m.faces().permute(permutation);
+}
+
+inline void optimize_vertices_for_faces(Mesh& m)
+{
+    uint64_t rng = 1;
+
+    std::vector<std::pair<int, int>> face_vertex_indices;
+    for (auto v : m.vertices())
+    {
+        face_handle ff;
+        auto cnt = 0;
+        for (auto f : v.faces())
+        {
+            if (f.is_invalid())
+                continue;
+
+            ++cnt;
+            if (detail::xorshift64star(rng) % cnt == 0)
+                ff = f;
+        }
+        face_vertex_indices.emplace_back(ff.idx.value, v.idx.value);
+    }
+
+    // sort by face idx
+    sort(face_vertex_indices.begin(), face_vertex_indices.end());
+
+    // extract vertex indices
+    std::vector<int> permutation(face_vertex_indices.size());
+    for (auto i = 0u; i < face_vertex_indices.size(); ++i)
+        permutation[i] = face_vertex_indices[i].second;
+
+    // apply permutation
+    m.vertices().permute(permutation);
+}
+}

src/polymesh/algorithms/properties.hh

@@ -101,16 +101,22 @@ Scalar angle_defect(vertex_handle v, vertex_attribute<Vec3> const& position);
 /// ======== IMPLEMENTATION ========
 
 inline bool is_boundary(vertex_handle v) { return v.is_boundary(); }
+inline bool is_vertex_boundary(vertex_handle v) { return v.is_boundary(); }
 
-inline bool is_boundary(face_handle v) { return v.is_boundary(); }
+inline bool is_boundary(face_handle f) { return f.is_boundary(); }
+inline bool is_face_boundary(face_handle f) { return f.is_boundary(); }
 
-inline bool is_boundary(edge_handle v) { return v.is_boundary(); }
+inline bool is_boundary(edge_handle e) { return e.is_boundary(); }
+inline bool is_edge_boundary(edge_handle e) { return e.is_boundary(); }
 
-inline bool is_boundary(halfedge_handle v) { return v.is_boundary(); }
+inline bool is_boundary(halfedge_handle h) { return h.is_boundary(); }
+inline bool is_halfedge_boundary(halfedge_handle h) { return h.is_boundary(); }
 
 inline bool is_isolated(vertex_handle v) { return v.is_isolated(); }
+inline bool is_vertex_isolated(vertex_handle v) { return v.is_isolated(); }
 
-inline bool is_isolated(edge_handle v) { return v.is_isolated(); }
+inline bool is_isolated(edge_handle e) { return e.is_isolated(); }
+inline bool is_edge_isolated(edge_handle e) { return e.is_isolated(); }
 
 inline int valence(vertex_handle v) { return v.adjacent_vertices().size(); }
 

src/polymesh/algorithms/stats.hh

@@ -6,6 +6,7 @@
 #include "../fields.hh"
 
 #include "components.hh"
+#include "properties.hh"
 
 namespace polymesh
 {
@@ -56,6 +57,8 @@ void print_stats(std::ostream& out, Mesh const& m, vertex_attribute<Vec3> const*
     // TODO: genus
     // TODO: boundaries
     // TODO: isolated verts, edges
+    out << "  Isolated Vertices: " << m.vertices().count(is_vertex_isolated);
+    out << "  Isolated Edges: " << m.edges().count(is_edge_isolated);
 
     if (position)
     {
@@ -72,6 +75,10 @@ void print_stats(std::ostream& out, Mesh const& m, vertex_attribute<Vec3> const*
 
         auto avg = m.vertices().avg(pos);
         out << "  Vertex Centroid: " << field_3d<Vec3>::to_string(avg) << ln;
+
+        auto el_minmax = m.edges().minmax([&](edge_handle e) { return edge_length(e, pos); });
+        auto el_avg = m.edges().avg([&](edge_handle e) { return edge_length(e, pos); });
+        out << "  Edge Lengths: " << el_minmax.min << " .. " << el_minmax.max << " (avg " << el_avg << ")" << ln;
     }
 }
 }

src/polymesh/attribute_base.hh

@@ -115,6 +115,7 @@ protected:
     primitive_attribute_base(Mesh const* mesh) : mMesh(mesh) {} // no registration, it's too early!
     virtual void on_resize(int new_size) = 0;
     virtual void apply_remapping(std::vector<int> const& map) = 0;
+    virtual void apply_transpositions(std::vector<std::pair<int, int>> const& ts) = 0;
     void register_attr();
     void deregister_attr();
     friend class Mesh;

src/polymesh/attributes.hh

@@ -76,6 +76,7 @@ protected:
 
     void on_resize(int newSize) override { mData.resize(newSize, mDefaultValue); }
     void apply_remapping(std::vector<int> const& map) override;
+    void apply_transpositions(std::vector<std::pair<int, int>> const& ts) override;
 
     // ctor
 protected:

src/polymesh/detail/permutation.hh

+#pragma once
+
+#include <valarray>
+#include <vector>
+
+namespace polymesh
+{
+namespace detail
+{
+/// Applies a permutation that is given by a remapping
+/// p[curr_idx] = new_idx
+/// Calculates the necessary transpositions and calls s(i, j) for each of it
+template <class Swap>
+void apply_permutation(std::vector<int> const& p, Swap&& s);
+
+/// Returns true if the parameter is actually a permutation
+bool is_valid_permutation(std::vector<int> const& p);
+
+/// Returns a list of transpositions that result in the given remapping
+/// p[curr_idx] = new_idx
+std::vector<std::pair<int, int>> transpositions_of(std::vector<int> const& p);
+
+/// ======== IMPLEMENTATION ========
+
+inline bool is_valid_permutation(std::vector<int> const& p)
+{
+    std::vector<int> r(p.size(), -1);
+
+    for (auto i = 0u; i < p.size(); ++i)
+    {
+        auto pi = p[i];
+        if (pi < 0 || pi >= p.size())
+            return false; // out of bound
+
+        if (r[pi] != -1)
+            return false; // not injective
+
+        r[pi] = i;
+    }
+
+    return true;
+}
+
+template <class Swap>
+void apply_permutation(std::vector<int> const& p, Swap&& s)
+{
+    auto size = p.size();
+    std::vector<int> cycle;
+    std::valarray<bool> visited(false, size);
+    for (auto pi = 0u; pi < size; ++pi)
+    {
+        auto i = pi;
+
+        if (visited[i])
+            continue;
+
+        cycle.clear();
+        visited[i] = true;
+        i = p[i];
+        while (!visited[i])
+        {
+            // mark
+            visited[i] = true;
+
+            // add to cycle
+            cycle.push_back(i);
+
+            // advance
+            i = p[i];
+        }
+
+        for (auto j = (int)cycle.size() - 1; j >= 0; --j)
+            s(pi, cycle[j]);
+    }
+}
+
+inline std::vector<std::pair<int, int>> transpositions_of(std::vector<int> const& p)
+{
+    std::vector<std::pair<int, int>> ts;
+    apply_permutation(p, [&](int i, int j) { ts.emplace_back(i, j); });
+    return ts;
+}
+}
+}

src/polymesh/detail/random.hh

+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+
+namespace polymesh
+{
+namespace detail
+{
+// from https://en.wikipedia.org/wiki/Xorshift
+uint64_t xorshift64star(uint64_t &state)
+{
+    uint64_t x = state; /* The state must be seeded with a nonzero value. */
+    x ^= x >> 12;       // a
+    x ^= x << 25;       // b
+    x ^= x >> 27;       // c
+    state = x;
+    return x * 0x2545F4914F6CDD1D;
+}
+}
+}

src/polymesh/impl/impl_attributes.hh

@@ -27,6 +27,14 @@ void primitive_attribute<tag, AttrT>::apply_remapping(const std::vector<int> &ma
         this->mData[i] = this->mData[map[i]];
 }
 
+template <class tag, class AttrT>
+void primitive_attribute<tag, AttrT>::apply_transpositions(std::vector<std::pair<int, int>> const &ts)
+{
+    using std::swap;
+    for (auto t : ts)
+        swap(this->mData[t.first], this->mData[t.second]);
+}
+
 template <class tag, class AttrT>
 primitive_attribute<tag, AttrT>::primitive_attribute(primitive_attribute const &rhs) noexcept : primitive_attribute_base<tag>(rhs.mMesh) // copy
 {

src/polymesh/impl/impl_mesh.hh

@@ -2,6 +2,8 @@
 
 #include "../Mesh.hh"
 
+#include "../detail/permutation.hh"
+
 namespace polymesh
 {
 inline vertex_index Mesh::add_vertex()
@@ -1020,6 +1022,99 @@ inline void Mesh::halfedge_rotate_prev(halfedge_index h)
     fix_boundary_state_of(h1_ref.face);
 }
 
+inline void Mesh::permute_vertices(std::vector<int> const &p)
+{
+    assert(detail::is_valid_permutation(p));
+
+    // calculate transpositions
+    auto ts = detail::transpositions_of(p);
+
+    // apply them
+    for (auto t : ts)
+        std::swap(mVertices[t.first], mVertices[t.second]);
+
+    // fix half-edges
+    for (auto &h : mHalfedges)
+        if (h.to_vertex.is_valid())
+            h.to_vertex.value = p[h.to_vertex.value];
+
+    // update attributes
+    for (auto a = mVertexAttrs; a; a = a->mNextAttribute)
+        a->apply_transpositions(ts);
+}
+
+inline void Mesh::permute_faces(std::vector<int> const &p)
+{
+    assert(detail::is_valid_permutation(p));
+
+    // calculate transpositions
+    auto ts = detail::transpositions_of(p);
+
+    // apply them
+    for (auto t : ts)
+        std::swap(mFaces[t.first], mFaces[t.second]);
+
+    // fix half-edges
+    for (auto &h : mHalfedges)
+        if (h.face.is_valid())
+            h.face.value = p[h.face.value];
+
+    // update attributes
+    for (auto a = mFaceAttrs; a; a = a->mNextAttribute)
+        a->apply_transpositions(ts);
+}
+
+inline void Mesh::permute_edges(std::vector<int> const &p)
+{
+    assert(detail::is_valid_permutation(p));
+
+    std::vector<int> hp(p.size() * 2);
+    for (auto i = 0u; i < p.size(); ++i)
+    {
+        hp[i * 2 + 0] = p[i] * 2 + 0;
+        hp[i * 2 + 1] = p[i] * 2 + 1;
+    }
+    assert(detail::is_valid_permutation(hp));
+
+    // calculate transpositions
+    std::vector<std::pair<int, int>> edge_ts;
+    std::vector<std::pair<int, int>> halfedge_ts;
+
+    detail::apply_permutation(p, [&](int i, int j) {
+        edge_ts.emplace_back(i, j);
+        halfedge_ts.emplace_back((i << 1) + 0, (j << 1) + 0);
+        halfedge_ts.emplace_back((i << 1) + 1, (j << 1) + 1);
+    });
+
+    // apply them
+    for (auto t : halfedge_ts)
+        std::swap(mHalfedges[t.first], mHalfedges[t.second]);
+
+    // fix half-edges
+    for (auto &v : mVertices)
+        if (v.outgoing_halfedge.value >= 0)
+            v.outgoing_halfedge.value = hp[v.outgoing_halfedge.value];
+
+    for (auto &f : mFaces)
+        if (f.halfedge.value >= 0)
+            f.halfedge.value = hp[f.halfedge.value];
+
+    for (auto &h : mHalfedges)
+    {
+        if (h.next_halfedge.value >= 0)
+            h.next_halfedge.value = hp[h.next_halfedge.value];
+
+        if (h.prev_halfedge.value >= 0)
+            h.prev_halfedge.value = hp[h.prev_halfedge.value];
+    }
+
+    // update attributes
+    for (auto a = mEdgeAttrs; a; a = a->mNextAttribute)
+        a->apply_transpositions(edge_ts);
+    for (auto a = mHalfedgeAttrs; a; a = a->mNextAttribute)
+        a->apply_transpositions(halfedge_ts);
+}
+
 inline void Mesh::compactify()
 {
     if (is_compact())

src/polymesh/impl/impl_ranges.hh

@@ -68,6 +68,17 @@ bool smart_range<this_t, ElementT>::all(PredT &&p) const
     return true;
 }
 
+template <class this_t, class ElementT>
+template <class PredT>
+int smart_range<this_t, ElementT>::count(PredT &&p) const
+{
+    auto cnt = 0;
+    for (auto h : *static_cast<this_t const *>(this))
+        if (p(h))
+            ++cnt;
+    return cnt;
+}
+
 template <class this_t, class ElementT>
 int smart_range<this_t, ElementT>::count() const
 {
@@ -435,6 +446,24 @@ void halfedge_collection<iterator>::remove_edge(halfedge_handle h) const
     this->mesh->remove_edge(this->mesh->edge_of(h.idx));
 }
 
+template <class iterator>
+void face_collection<iterator>::permute(std::vector<int> const &p) const
+{
+    this->mesh->permute_faces(p);
+}
+
+template <class iterator>
+void edge_collection<iterator>::permute(std::vector<int> const &p) const
+{
+    this->mesh->permute_edges(p);
+}
+
+template <class iterator>
+void vertex_collection<iterator>::permute(std::vector<int> const &p) const
+{
+    this->mesh->permute_vertices(p);
+}
+
 
 inline valid_vertex_collection Mesh::vertices()
 {

src/polymesh/ranges.hh

@@ -51,6 +51,9 @@ struct smart_range
     /// NOTE: this is an O(n) operation, prefer size() if available
     /// TODO: maybe SFINAE to implement this via size() if available?
     int count() const;
+    /// returns the number of elements fulfilling p(v) in this range
+    template <class PredT>
+    int count(PredT&& p) const;
 
     /// calculates min(f(e)) over all elements
     /// undefined behavior if range is empty
@@ -184,6 +187,11 @@ struct vertex_collection : smart_collection<Mesh*, vertex_tag, iterator>
     /// Removes a vertex (and all adjacent faces and edges)
     /// (marks them as removed, compactify mesh to actually remove them)
     void remove(vertex_handle v) const;
+
+    /// applies an index remapping to all vertex indices
+    /// p[curr_idx] = new_idx
+    /// NOTE: invalidates all affected handles/iterators
+    void permute(std::vector<int> const& p) const;
 };
 
 /// Collection of all faces of a mesh
@@ -219,6 +227,11 @@ struct face_collection : smart_collection<Mesh*, face_tag, iterator>
     /// Removes a face (adjacent edges and vertices are NOT removed)
     /// (marks it as removed, compactify mesh to actually remove it)
     void remove(face_handle f) const;
+
+    /// applies an index remapping to all face indices
+    /// p[curr_idx] = new_idx
+    /// NOTE: invalidates all affected handles/iterators
+    void permute(std::vector<int> const& p) const;
 };
 
 /// Collection of all edges of a mesh
@@ -251,6 +264,11 @@ struct edge_collection : smart_collection<Mesh*, edge_tag, iterator>
     /// Removes an edge (and both adjacent faces, vertices are NOT removed)
     /// (marks them as removed, compactify mesh to actually remove them)
     void remove(edge_handle e) const;
+
+    /// applies an index remapping to all edge indices