toddcox-visualize/vis/include/solver.hpp

#pragma once

#include <tc/core.hpp>
#include <cmath>
#include <optional>
#include <numeric>
#include <iostream>

#include <geometry.hpp>

/**
 * 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<int N>
[[nodiscard]]
Indices<N> recontext(
    Indices<N> prims,
    const tc::Group &G,
    const tc::Group &H
) {
    const auto table = G.solve({});
    const auto cosets = H.solve({});

    tc::Path<size_t> path(cosets, H.gens());

    std::vector<size_t> map(path.order());
    path.walk(0, [&table](size_t coset, size_t gen) {
        return table.get(coset, gen);
    }, map.begin());

    Indices<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<int N>
Indices<N> merge(
    const std::vector<Indices<N>> &meshes
) {
    size_t cols = 0;
    for (const auto &mesh: meshes) {
        cols += mesh.cols();
    }

    Indices<N> res(N, cols);

    size_t offset = 0;
    for (const Indices<N> &mesh: meshes) {
        res.middleCols(offset, mesh.cols()) = mesh;
        offset += mesh.cols();
    }

    return res;
}

template<int N>
[[nodiscard]]
std::vector<Indices<N>> tile(
    Indices<N> prims,
    const tc::Group &G,
    const tc::Group &H
) {
    const auto &table = G.solve({});
    const auto &cosets = G.solve(H.gens());

    tc::Path<> path(cosets);

    std::vector<Indices<N>> res(path.order());

    path.walk(prims, [&](Indices<N> from, auto gen) {
        for (int i = 0; i < from.size(); ++i) {
            from(i) = table.get(from(i), gen);
        }
        return from;
    }, res.begin());

    return res;
}

/**
 * Produce a mesh of higher dimension by fanning a single point to all primitives in this mesh.
 */
template<int N>
[[nodiscard]]
Indices<N + 1> fan(
    Indices<N> prims,
    int root
) {
    Indices<N + 1> res(N + 1, prims.cols());

    res.topRows(N) = prims;
    res.bottomRows(1) = root;

    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<int N>
Indices<N> cell(
    const tc::Group &G
) {
    assert(G.rank() + 1 == N);  // g_gens size must be one less than N

    std::vector<Indices<N - 1>> facets;

    for (auto H: G.subs(N - 2)) {
        auto sub_base = cell<N - 1>(H);
        auto base = recontext(sub_base, G, H);
        auto tiles = tile(base, G, H);

        facets.insert(
            facets.end(),
            tiles.begin() + 1, // skip "near" facet
            tiles.end()
        );
    }

    return fan(merge(facets), 0);
}

/**
 * Single-index primitives should not be further triangulated.
 */
template<>
Indices<1> cell<1>(
    const tc::Group &G
) {
    assert(G.rank() == 0); // rank must be 0 for trivial Mesh

    return Indices<1>::Zero(1, 1);
}

template<int N>
struct Hull {
    struct Tiling {
        Eigen::Index first;
        Eigen::Index count;
    };

    std::vector<Tiling> tilings{};
    std::vector<std::vector<size_t>> subgroups{};

    Indices<N> inds{};

    explicit Hull(tc::Group const &G) {
        std::vector<Indices<N>> parts;

        Eigen::Index first = 0;
        for (const auto &H: G.subs(N - 1)) {
            auto sub_base = cell<N>(H);
            auto base = recontext(sub_base, G, H);
            auto tiles = tile(base, G, H);
            auto tiling = merge(tiles);

            subgroups.push_back(H.gens());
            tilings.push_back({first, tiling.cols()});

            parts.push_back(tiling);
            first += tiling.cols();
        }

        inds = merge(parts);
    }
};

struct Points {
    Eigen::Array<float, 4, Eigen::Dynamic> verts;

    explicit Points(tc::Group const &G, Eigen::Vector<float, 5> root) {
        auto mirrors = mirror<5>(G);
        auto corners = plane_intersections(mirrors);

        auto start = corners * root;

        tc::Cosets table = G.solve();
        tc::Path<vec5> path(table, mirrors.colwise());

        Eigen::Array<float, 5, Eigen::Dynamic> higher(5, path.order());
        path.walk(start, Reflect(), higher.matrix().colwise().begin());

        verts = Stereo()(higher);
    }
};