WIP: Move old geometry solver to include/geo

include/geo/combo.hpp

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <set>
+#include <algorithm>
+
+template<typename V, typename M>
+V select(const V &data, const M &mask, size_t count) {
+    V result;
+    result.reserve(count);
+
+    for (int i = 0; i < mask.size(); ++i) {
+        if (mask[i]) result.push_back(data[i]);
+    }
+
+    return result;
+}
+
+template<typename V>
+std::vector<V> combinations(const V &data, const size_t count) {
+    std::vector<V> result;
+
+    std::vector<bool> mask(data.size(), false);
+    std::fill(mask.begin(), mask.begin() + count, true);
+
+    do {
+        result.push_back(select(data, mask, count));
+    } while (std::next_permutation(mask.begin(), mask.end(), std::greater<>()));
+
+    return result;
+}

include/geo/geometry.hpp

@@ -0,0 +1,63 @@
+#pragma once
+
+#include <tc/core.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+
+#include <Eigen/Eigen>
+
+#include "combo.hpp"
+
+template<unsigned N>
+using Prims = Eigen::Matrix<unsigned, N, Eigen::Dynamic>;
+
+template<int N>
+using vec = Eigen::Matrix<float, N, 1>;
+template<int N>
+using mat = Eigen::Matrix<float, N, N>;
+
+using vec1 = vec<1>;
+using vec2 = vec<2>;
+using vec3 = vec<3>;
+using vec4 = vec<4>;
+using vec5 = vec<5>;
+
+using mat1 = mat<1>;
+using mat2 = mat<2>;
+using mat3 = mat<3>;
+using mat4 = mat<4>;
+using mat5 = mat<5>;
+
+mat4 orthographic(float left, float right, float bottom, float top, float front, float back) {
+    mat4 res = mat4();
+    res <<
+        2 / (right - left), 0, 0, -(right + left) / (right - left),
+            0, 2 / (top - bottom), 0, -(top + bottom) / (top - bottom),
+            0, 0, 2 / (front - back), -(front + back) / (front - back),
+            0, 0, 0, 1;
+    return res;
+}
+
+mat4 perspective(float fovy, float aspect, float zNear, float zFar) {
+    float tanHalfFovy(std::tan(fovy / 2));
+
+    mat4 res = mat4::Identity();
+    res(0, 0) = 1 / (aspect * tanHalfFovy);
+    res(1, 1) = 1 / (tanHalfFovy);
+    res(2, 2) = -(zFar + zNear) / (zFar - zNear);
+    res(3, 2) = -1;
+    res(2, 3) = -(2 + zFar * zNear) / (zFar - zNear);
+    return res;
+}
+
+mat4 translation(float x, float y, float z) {
+    mat4 res = mat4();
+    res <<
+        1, 0, 0, x,
+        0, 1, 0, y,
+        0, 0, 1, z,
+        0, 0, 0, 1;
+    return res;
+}

include/geo/mirror.hpp

@@ -0,0 +1,128 @@
+#pragma once
+
+#include <tc/core.hpp>
+
+#include <cmath>
+#include <vector>
+#include <algorithm>
+
+#include <geo/geometry.hpp>
+
+template<class V>
+float dot(int n, const V &a, const V &b) {
+    float sum = 0;
+    for (int i = 0; i < n; ++i) {
+        sum += a[i] * b[i];
+    }
+    return sum;
+}
+
+template<unsigned N>
+std::vector<vec<N>> mirror(const tc::Group &group) {
+    std::vector<std::vector<float>> mirrors;
+
+    for (int p = 0; p < group.ngens; ++p) {
+        std::vector<float> vp;
+        for (int m = 0; m < p; ++m) {
+            auto &vq = mirrors[m];
+            vp.push_back((cos(M_PI / group.get(p, m)) - dot(m, vp, vq)) / vq[m]);
+        }
+        vp.push_back(std::sqrt(1 - dot(p, vp, vp)));
+
+        for (const auto &v : mirrors) {
+            if (dot(p, vp, vp) > 0) {
+                for (auto &e : vp) {
+                    e *= -1;
+                }
+                break;
+            }
+        }
+
+        mirrors.push_back(vp);
+    }
+
+    std::vector<vec<N>> res;
+    for (const auto &v : mirrors) {
+        vec<N> rv = vec<N>::Zero();
+
+        // ortho proj
+        for (int i = 0; i < std::min(v.size(), (size_t) N); ++i) {
+            rv[i] = v[i];
+        }
+
+        res.push_back(rv);
+    }
+    return res;
+}
+
+template<unsigned N>
+vec<N> stereo(const vec<N + 1> &v) {
+    vec<N> r;
+    for (int i = 0; i < N; ++i) {
+        r[i] = v[i] / (1 - v[N]);
+    }
+    return r;
+}
+
+template<unsigned N>
+vec<N> ortho(const vec<N + 1> &v) {
+    vec<N> r;
+    for (int i = 0; i < N; ++i) {
+        r[i] = v[i];
+    }
+    return r;
+}
+
+template<class V>
+V project(const V &vec, const V &target) {
+    return vec.dot(target) / target.dot(target) * target;
+}
+
+template<class V>
+V reflect(const V &a, const V &axis) {
+    return a - 2.f * project(a, axis);
+}
+
+template<class V>
+V gram_schmidt_last(std::vector<V> vecs) {
+    for (int i = 0; i < vecs.size(); ++i) {
+        for (int j = 0; j < i; ++j) {
+            vecs[i] -= project(vecs[i], vecs[j]);
+        }
+    }
+
+    return vecs[vecs.size() - 1].normalized();
+}
+
+template<class V, class C>
+V barycentric(const std::vector<V> &basis, const C &coords) {
+    V res = V::Zero();
+
+    int N = std::min((int) basis.size(), (int) coords.rows());
+    for (int i = 0; i < N; ++i) {
+        res += basis[i] * coords[i];
+    }
+    return res;
+}
+
+template<class V>
+std::vector<V> plane_intersections(std::vector<V> normals) {
+    std::vector<V> results(normals.size());
+
+    for (int i = 0; i < normals.size(); ++i) {
+        std::rotate(normals.begin(), normals.begin() + 1, normals.end());
+        results[i] = gram_schmidt_last(normals);
+    }
+
+    return results;
+}
+
+template<unsigned N>
+mat<N> rot(int u, int v, float theta) {
+    mat<N> res = mat<N>::Identity();
+    res(u, u) = std::cos(theta);
+    res(u, v) = std::sin(theta);
+    res(v, u) = -std::sin(theta);
+    res(v, v) = std::cos(theta);
+    return res;
+}

include/geo/solver.hpp

@@ -0,0 +1,233 @@
+#pragma once
+
+#include <tc/core.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+
+#include <geo/geometry.hpp>
+
+#include "combo.hpp"
+
+/**
+ * Produce a list of all generators for the group context. The range [0..group.ngens).
+ */
+std::vector<int> generators(const tc::Group &context) {
+    // todo if tc::Group has 'global' generators, then this will be a member of tc::Group.
+    //  std::iota would populate a 'default' list of names, if names are not provided.
+    std::vector<int> g_gens(context.ngens);
+    std::iota(g_gens.begin(), g_gens.end(), 0);
+    return g_gens;
+}
+
+/**
+ * Determine which of g_gens are the correct names for sg_gens within the current context
+ */
+std::vector<int> recontext_gens(
+    const tc::Group &context,
+    std::vector<int> g_gens,
+    std::vector<int> sg_gens) {
+    // todo ideally tc::Group will deal in 'global' generators so this stell will be unecessary.
+
+    std::sort(g_gens.begin(), g_gens.end());
+
+    int inv_gen_map[context.ngens];
+    for (size_t i = 0; i < g_gens.size(); i++) {
+        inv_gen_map[g_gens[i]] = i;
+    }
+
+    std::vector<int> s_sg_gens;
+    s_sg_gens.reserve(sg_gens.size());
+    for (const auto gen: sg_gens) {
+        s_sg_gens.push_back(inv_gen_map[gen]);
+    }
+    std::sort(s_sg_gens.begin(), s_sg_gens.end());
+
+    return s_sg_gens;
+}
+
+/**
+ * Solve the cosets generated by sg_gens within the subgroup generated by g_gens of the group context
+ */
+tc::Cosets solve(
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    // todo this should also be handled with 'global' generators.
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+    return context.subgroup(g_gens).solve(proper_sg_gens);
+}
+
+/**
+ * Apply some context transformation to all primitives of this mesh.
+ */
+template<unsigned N>
+void apply(const tc::Cosets &table, int gen, Prims<N> &mat) {
+    auto data = mat.data();
+    for (int i = 0; i < mat.size(); ++i) {
+        data[i] = table.get(data[i], gen);
+    }
+}
+
+/**
+ * Convert the indexes of this mesh to those of a different context, using g_gens to build the parent context and sg_gens to build this context.
+ */
+template<unsigned N>
+[[nodiscard]]
+Prims<N> recontext(
+    Prims<N> prims,
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    // todo this will be simpler with 'global' gens, but it's still not free...
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+    const auto table = solve(context, g_gens, {});
+    const auto path = solve(context, sg_gens, {}).path;
+
+    auto map = path.template walk<int, int>(0, proper_sg_gens, [table](int coset, int gen) {
+        return table.get(coset, gen);
+    });
+
+    Prims<N> res(prims);
+    auto data = res.data();
+    for (int i = 0; i < prims.size(); ++i) {
+        data[i] = map[data[i]];
+    }
+
+    return res;
+}
+
+/**
+ * Union several meshes of the same dimension
+ */
+template<unsigned N>
+Prims<N> merge(const std::vector<Prims<N>> &meshes) {
+    // todo (?) might be possible with NullaryExpr
+    size_t cols = 0;
+    for (const auto &mesh: meshes) {
+        cols += mesh.cols();
+    }
+
+    Prims<N> res(N, cols);
+
+    size_t offset = 0;
+    for (const Prims<N> &mesh: meshes) {
+        res.middleCols(offset, mesh.cols()) = mesh;
+        offset += mesh.cols();
+    }
+
+    return res;
+}
+
+template<unsigned N>
+[[nodiscard]]
+std::vector<Prims<N>> tile(
+    Prims<N> prims,
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    // todo convert to nullaryexpr.
+    //  some stuff will be easier with global generators, but not all.
+    Prims<N> base = recontext<N>(prims, context, g_gens, sg_gens);
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+
+    const auto table = solve(context, g_gens, {});
+    const auto path = solve(context, g_gens, sg_gens).path;
+
+    std::vector<int> _gens = generators(context);
+
+    std::vector<Prims<N>> res = path.walk<Prims<N>, int>(
+        base, _gens,
+        [&](Prims<N> from, int gen) {
+            apply<N>(table, gen, from);
+            return from;
+        }
+    );
+
+    return res;
+}
+
+/**
+ * Produce a mesh of higher dimension by fanning a single point to all primitives in this mesh.
+ */
+template<unsigned N>
+[[nodiscard]]
+Prims<N + 1> fan(Prims<N> prims, int root) {
+    // todo convert to nullaryexpr.
+    Prims<N + 1> res(N + 1, prims.cols());
+
+    res.topRows(1) = Prims<1>::Constant(1, prims.cols(), root);
+    res.bottomRows(N) = prims;
+
+    return res;
+}
+
+/**
+ * Produce a mesh of primitives that fill out the volume of the subgroup generated by generators g_gens within the group context
+ */
+template<unsigned N>
+Prims<N> triangulate(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    // todo (?) might be possible with nullaryexpr
+    //  not so sure, though.
+    if (g_gens.size() + 1 != N) // todo make static assert
+        throw std::logic_error("g_gens size must be one less than N");
+
+    const auto &combos = combinations(g_gens, g_gens.size() - 1);
+
+    std::vector<Prims<N>> meshes;
+
+    for (const auto &sg_gens: combos) {
+        auto base = triangulate<N - 1>(context, sg_gens);
+        auto parts = tile<N - 1>(base, context, g_gens, sg_gens);
+        parts.erase(parts.begin(), parts.begin() + 1);
+        auto raised = merge<N - 1>(parts);
+        auto fanned = fan<N - 1>(raised, 0);
+        meshes.push_back(fanned);
+    }
+
+    return merge<N>(meshes);
+}
+
+/**
+ * Single-index primitives should not be further triangulated.
+ */
+template<>
+Prims<1> triangulate<1>(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    if (not g_gens.empty()) // todo make static assert
+        throw std::logic_error("g_gens must be empty for a trivial Mesh");
+
+    return Prims<1>::Zero(1, 1);
+}
+
+template<unsigned N, class T>
+auto hull(const tc::Group &group, T all_sg_gens, const std::vector<std::vector<int>> &exclude) {
+    std::vector<Prims<N>> parts;
+    auto g_gens = generators(group);
+    for (const std::vector<int> &sg_gens: all_sg_gens) {
+        bool excluded = false;
+        for (const auto &test: exclude) {
+            if (sg_gens == test) {
+                excluded = true;
+                break;
+            }
+        }
+        if (excluded) continue;
+
+        const auto &base = triangulate<N>(group, sg_gens);
+        const auto &tiles = tile<N>(base, group, g_gens, sg_gens);
+        for (const auto &tile: tiles) {
+            parts.push_back(tile);
+        }
+    }
+    return parts;
+}