Multiple passes for different point-cloud resolutions

.gitignore

@@ -1,152 +0,0 @@
-# Created by .ignore support plugin (hsz.mobi)
-### C template
-# Prerequisites
-*.d
-
-# Object files
-*.o
-*.ko
-*.obj
-*.elf
-
-# Linker output
-*.ilk
-*.map
-*.exp
-
-# Precompiled Headers
-*.gch
-*.pch
-
-# Libraries
-*.lib
-*.a
-*.la
-*.lo
-
-# Shared objects (inc. Windows DLLs)
-*.dll
-*.so
-*.so.*
-*.dylib
-
-# Executables
-*.exe
-*.out
-*.app
-*.i*86
-*.x86_64
-*.hex
-
-# Debug files
-*.dSYM/
-*.su
-*.idb
-*.pdb
-
-# Kernel Module Compile Results
-*.mod*
-*.cmd
-.tmp_versions/
-modules.order
-Module.symvers
-Mkfile.old
-dkms.conf
-### CMake template
-CMakeCache.txt
-CMakeFiles
-CMakeScripts
-Testing
-Makefile
-cmake_install.cmake
-install_manifest.txt
-compile_commands.json
-CTestTestfile.cmake
-### JetBrains template
-# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio and WebStorm
-# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839
-
-# User-specific stuff
-.idea/**/workspace.xml
-.idea/**/tasks.xml
-.idea/**/dictionaries
-.idea/**/shelf
-
-# Sensitive or high-churn files
-.idea/**/dataSources/
-.idea/**/dataSources.ids
-.idea/**/dataSources.local.xml
-.idea/**/sqlDataSources.xml
-.idea/**/dynamic.xml
-.idea/**/uiDesigner.xml
-.idea/**/dbnavigator.xml
-
-# Gradle
-.idea/**/gradle.xml
-.idea/**/libraries
-
-# CMake
-cmake-build-debug/
-cmake-build-release/
-
-# Mongo Explorer plugin
-.idea/**/mongoSettings.xml
-
-# File-based project format
-*.iws
-
-# IntelliJ
-out/
-
-# mpeltonen/sbt-idea plugin
-.idea_modules/
-
-# JIRA plugin
-atlassian-ide-plugin.xml
-
-# Cursive Clojure plugin
-.idea/replstate.xml
-
-# Crashlytics plugin (for Android Studio and IntelliJ)
-com_crashlytics_export_strings.xml
-crashlytics.properties
-crashlytics-build.properties
-fabric.properties
-
-# Editor-based Rest Client
-.idea/httpRequests
-### C++ template
-# Prerequisites
-*.d
-
-# Compiled Object files
-*.slo
-*.lo
-*.o
-*.obj
-
-# Precompiled Headers
-*.gch
-*.pch
-
-# Compiled Dynamic libraries
-*.so
-*.dylib
-*.dll
-
-# Fortran module files
-*.mod
-*.smod
-
-# Compiled Static libraries
-*.lai
-*.la
-*.a
-*.lib
-
-# Executables
-*.exe
-*.out
-*.app
-
-.idea/

.gitmodules

@@ -1,7 +0,0 @@
-[submodule "vendor/glfw"]
-	path = vendor/glfw
-	url = https://github.com/glfw/glfw.git
-[submodule "vendor/glad"]
-	path = vendor/glad
-	url = https://github.com/Dav1dde/glad.git
-	branch = c

CMakeLists.txt

@@ -1,14 +1,31 @@
-cmake_minimum_required(VERSION 3.10)
+cmake_minimum_required(VERSION 3.21)
 project(hopf)
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 
-add_library(glad vendor/glad/src/glad.c)
-target_include_directories(glad PUBLIC vendor/glad/include)
+include(FetchContent)
 
-option(GLFW_BUILD_DOCS OFF)
-option(GLFW_BUILD_EXAMPLES OFF)
-option(GLFW_BUILD_TESTS OFF)
-add_subdirectory(vendor/glfw)
+include(External/glfw.cmake)
+include(External/imgui.cmake)
+include(External/eigen.cmake)
+include(External/glad.cmake)
+include(External/json.cmake)
 
-add_subdirectory(main)
+include_directories(include)
+
+add_custom_target(resources DEPENDS resources_output)
+add_custom_command(
+    OUTPUT resources_output
+    COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/res ${CMAKE_CURRENT_BINARY_DIR}/res
+    COMMENT "Copying Resources")
+
+add_executable(hopf
+    src/main.cpp
+    src/gl/debug.hpp
+    src/gl/shader.hpp
+    src/gl/buffer.hpp
+    src/gl/vertexarray.hpp
+    src/gl/types.hpp)
+
+target_link_libraries(hopf glfw glad imgui eigen nlohmann_json)
+add_dependencies(hopf resources)

External/eigen.cmake

@@ -0,0 +1,9 @@
+FetchContent_Declare(
+    eigen
+    GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+    GIT_TAG 3.4
+)
+set(EIGEN_BUILD_DOC OFF CACHE INTERNAL "")
+set(BUILD_TESTING OFF CACHE INTERNAL "")
+set(EIGEN_BUILD_PKGCONFIG OFF CACHE INTERNAL "")
+FetchContent_MakeAvailable(eigen)

External/glad.cmake

@@ -0,0 +1,9 @@
+FetchContent_Declare(
+    glad
+    GIT_REPOSITORY https://github.com/Dav1dde/glad.git
+    GIT_TAG v0.1.36
+)
+set(GLAD_PROFILE "core" CACHE INTERNAL "OpenGL profile")
+set(GLAD_API "gl=4.6" CACHE INTERNAL "API type/version pairs, like \"gl=3.2,gles=\", no version means latest")
+set(GLAD_GENERATOR "c" CACHE INTERNAL "Language to generate the binding for")
+FetchContent_MakeAvailable(glad)

External/glfw.cmake

@@ -0,0 +1,10 @@
+FetchContent_Declare(
+    glfw
+    GIT_REPOSITORY https://github.com/glfw/glfw.git
+    GIT_TAG 3.3.6
+)
+set(GLFW_BUILD_EXAMPLES OFF CACHE INTERNAL "")
+set(GLFW_BUILD_TESTS OFF CACHE INTERNAL "")
+set(GLFW_BUILD_DOCS OFF CACHE INTERNAL "")
+set(GLFW_INSTALL OFF CACHE INTERNAL "")
+FetchContent_MakeAvailable(glfw)

External/imgui.cmake

@@ -0,0 +1,22 @@
+FetchContent_Declare(
+    imgui
+    GIT_REPOSITORY https://github.com/ocornut/imgui
+    GIT_TAG v1.86
+)
+FetchContent_MakeAvailable(imgui)
+
+add_library(
+    imgui
+    ${imgui_SOURCE_DIR}/backends/imgui_impl_glfw.cpp
+    ${imgui_SOURCE_DIR}/backends/imgui_impl_glfw.h
+    ${imgui_SOURCE_DIR}/backends/imgui_impl_opengl3.cpp
+    ${imgui_SOURCE_DIR}/backends/imgui_impl_opengl3.h
+    ${imgui_SOURCE_DIR}/imgui.cpp
+    ${imgui_SOURCE_DIR}/imgui.h
+    ${imgui_SOURCE_DIR}/imgui_demo.cpp
+    ${imgui_SOURCE_DIR}/imgui_draw.cpp
+    ${imgui_SOURCE_DIR}/imgui_tables.cpp
+    ${imgui_SOURCE_DIR}/imgui_widgets.cpp
+)
+target_include_directories(imgui PUBLIC ${imgui_SOURCE_DIR})
+target_link_libraries(imgui PRIVATE glfw)

External/json.cmake

@@ -0,0 +1,7 @@
+FetchContent_Declare(
+    json
+    GIT_REPOSITORY https://github.com/nlohmann/json.git
+    GIT_TAG v3.10.5
+)
+SET(JSON_ImplicitConversions OFF CACHE INTERNAL "")
+FetchContent_MakeAvailable(json)

README.md

@@ -1,21 +0,0 @@
-# hopf-fibration
-
-Generate images of the hopf fibration.
-
-![zoomed](images/zoomed.png)
-
-## Usage
-
-Hold keys `Q`, `W`, `E`, `A`, `S`, `D` to select from the six planes of rotation, and use the scroll wheel to rotate the fibration through the selected plane. The resulting Hopf links are stereographically projected into 3D space, and the resulting circles are rendered orthographically on screen. 
-
-Hold left-control `LCTL` and scroll to zoom in on the fibration. This affects only the orthographic rendering, not the fibration itself.
-
-![default](images/default.png)
-
-![rotated](images/rotated.png)
-
----
-
-## More information
-
-[Wikipedia: Hopf Fibration](https://en.wikipedia.org/wiki/Hopf_fibration)

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;
+}

include/ml/meshlib.hpp

@@ -0,0 +1,83 @@
+#pragma once
+
+#include <iostream>
+#include <ostream>
+#include <memory>
+#include <utility>
+
+#include <Eigen/Eigen>
+
+namespace ml {
+    using Matrix1Xui = Eigen::Matrix<unsigned int, 1, Eigen::Dynamic>;
+    using Matrix2Xui = Eigen::Matrix<unsigned int, 2, Eigen::Dynamic>;
+    using Matrix3Xui = Eigen::Matrix<unsigned int, 3, Eigen::Dynamic>;
+    using Matrix4Xui = Eigen::Matrix<unsigned int, 4, Eigen::Dynamic>;
+
+    template<typename PT_, typename CT_>
+    class Mesh {
+    public:
+        using Points = PT_;
+        using Cells = CT_;
+
+        Points points;
+        Cells cells;
+
+        Mesh(Points points, Cells cells)
+            : points(std::move(points)), cells(std::move(cells)) {}
+    };
+
+    auto make_cube(float radius) {
+        Eigen::Matrix3Xf points(3, 8);
+        points.fill(radius);
+        for (int i = 0; i < points.cols(); ++i) {
+            for (int j = 0; j < 3; ++j) {
+                if ((i >> j) & 1) {
+                    points(j, i) *= -1;
+                }
+            }
+        }
+
+        Matrix3Xui cells(3, 12);
+        cells.transpose()
+            << 0b000, 0b001, 0b010, 0b001, 0b010, 0b011,
+            0b100, 0b101, 0b110, 0b101, 0b110, 0b111,
+
+            0b000, 0b001, 0b100, 0b001, 0b100, 0b101,
+            0b010, 0b011, 0b110, 0b011, 0b110, 0b111,
+
+            0b000, 0b010, 0b100, 0b010, 0b100, 0b110,
+            0b001, 0b011, 0b101, 0b011, 0b101, 0b111;
+
+        return Mesh(points, cells);
+    }
+    
+    template<size_t Dim>
+    auto make_cube_wire(float radius) {
+        constexpr size_t NPoints = 1 << Dim;
+        constexpr size_t NCells = Dim * (NPoints >> 1);
+
+        Eigen::Matrix<float, Dim, NPoints> points;
+        points.fill(radius);
+        for (int i = 0; i < points.cols(); ++i) {
+            for (int j = 0; j < Dim; ++j) {
+                if ((i >> j) & 1) {
+                    points(j, i) *= -1;
+                }
+            }
+        }
+
+        Eigen::Matrix<unsigned int, 2, NCells> cells;
+        int k = 0;
+        for (int i = 0; i < NPoints; ++i) {
+            for (int j = 0; j < Dim; ++j) {
+                if ((i >> j) & 1) {
+                    cells(0, k) = i;
+                    cells(1, k) = i ^ (1 << j);
+                    k++;
+                }
+            }
+        }
+
+        return Mesh(points, cells);
+    }
+}

include/ml/meshlib_json.hpp

@@ -0,0 +1,68 @@
+#pragma once
+
+#include "meshlib.hpp"
+
+#include <Eigen/Eigen>
+#include <nlohmann/json.hpp>
+
+namespace Eigen {
+    template<class Derived>
+    void to_json(nlohmann::json &json, const Eigen::PlainObjectBase<Derived> &mat) {
+        using Scalar = typename Derived::Scalar;
+
+        auto rows = mat.rows();
+        auto cols = mat.cols();
+
+        std::vector<Scalar> vals(mat.size());
+        Map<Derived>(vals.data(), rows, cols) = mat;
+
+        json = {
+            {"rows", rows},
+            {"cols", cols},
+            {"vals", vals},
+        };
+    }
+
+    template<class Derived>
+    void from_json(const nlohmann::json &j, Derived &d) {
+        using Scalar = typename Derived::Scalar;
+
+        auto rows = j["rows"].get<Index>();
+        auto cols = j["cols"].get<Index>();
+        auto vals = j["vals"].get<std::vector<Scalar>>();
+
+        d = Map<Derived>(vals.data(), rows, cols);
+    }
+}
+
+namespace nlohmann {
+    template<typename PT_, typename CT_>
+    struct adl_serializer<ml::Mesh<PT_, CT_>> {
+        static void to_json(json &j, const ml::Mesh<PT_, CT_> &m) {
+            j = {
+                {"points", m.points},
+                {"cells",  m.cells},
+            };
+        }
+
+        static ml::Mesh<PT_, CT_> from_json(const json &j) {
+            return ml::Mesh<PT_, CT_>(
+                j["points"].get<PT_>(),
+                j["cells"].get<CT_>()
+            );
+        }
+    };
+}
+
+namespace ml {
+    template<typename PT_, typename CT_>
+    void write(const ml::Mesh<PT_, CT_> &mesh, std::ostream &&out) {
+        nlohmann::json json = mesh;
+        nlohmann::json::to_msgpack(json, out);
+    }
+
+    template<typename M_>
+    M_ read(std::istream &&in) {
+        return nlohmann::json::from_msgpack(in).get<M_>();
+    }
+}

main/CMakeLists.txt

@@ -1,29 +0,0 @@
-project(main)
-
-add_executable(${PROJECT_NAME}
-        src/main.cpp)
-
-target_link_libraries(${PROJECT_NAME}
-        glad
-        glfw)
-
-target_include_directories(${PROJECT_NAME}
-        PRIVATE include)
-
-set(SHADERS
-        shaders/main.frag
-        shaders/main.vert)
-add_custom_target(shaders DEPENDS ${SHADERS})
-
-add_custom_command(
-        TARGET shaders PRE_BUILD
-        COMMAND ${CMAKE_COMMAND} -E remove_directory ${CMAKE_CURRENT_BINARY_DIR}/shaders/
-        COMMENT "clearing shaders"
-)
-add_custom_command(
-        TARGET shaders POST_BUILD
-        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/shaders/ ${CMAKE_CURRENT_BINARY_DIR}/shaders/
-        COMMENT "copying shaders"
-)
-
-add_dependencies(${PROJECT_NAME} shaders)

main/include/glga.hpp

@@ -1,165 +0,0 @@
-#include <math.h>
-
-namespace ga {
-    struct Vec;
-    struct Bivec;
-    struct Mat;
-
-    struct Vec {
-        float x, y, z, w;
-
-        Vec(float x, float y, float z, float w) :
-            x(x), y(y), z(z), w(w) {}
-    };
-
-    struct Bivec {
-        float xy, yz, zw, wx, xz, yw;
-
-        Bivec(float xy, float yz, float zw, float wx, float xz, float yw) :
-            xy(xy), yz(yz), zw(zw), wx(wx), xz(xz), yw(yw) {}
-    };
-
-    struct Mat {
-        Vec x, y, z, w;
-
-        Mat(Vec x, Vec y, Vec z, Vec w) : x(x), y(y), z(z), w(w) {}
-    };
-
-    namespace unit {
-        Vec x() { return Vec(1, 0, 0, 0); };
-
-        Vec y() { return Vec(0, 1, 0, 0); };
-
-        Vec z() { return Vec(0, 0, 1, 0); };
-
-        Vec w() { return Vec(0, 0, 0, 1); };
-
-        Bivec xy() { return Bivec(1, 0, 0, 0, 0, 0); }
-
-        Bivec yz() { return Bivec(0, 1, 0, 0, 0, 0); }
-
-        Bivec zw() { return Bivec(0, 0, 1, 0, 0, 0); }
-
-        Bivec wx() { return Bivec(0, 0, 0, 1, 0, 0); }
-
-        Bivec xz() { return Bivec(0, 0, 0, 0, 1, 0); }
-
-        Bivec yw() { return Bivec(0, 0, 0, 0, 0, 1); }
-
-        Mat identity() { return Mat(x(), y(), z(), w()); }
-    }
-
-    float length2(Vec v) {
-        return v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w;
-    }
-
-    float length2(Bivec v) {
-        return v.xy * v.xy + v.yz * v.yz + v.zw * v.zw
-            + v.wx * v.wx + v.xz * v.xz + v.yw * v.yw;
-    }
-
-    float length(Vec v) { return sqrt(length2(v)); }
-
-    float length(Bivec v) { return sqrt(length2(v)); }
-
-    Vec add(Vec u, Vec v) {
-        return Vec(u.x + v.x,
-            u.y + v.y,
-            u.z + v.z,
-            u.w + v.w);
-    }
-
-    Bivec add(Bivec u, Bivec v) {
-        Bivec(u.xy + v.xy,
-            u.yz + v.yz,
-            u.zw + v.zw,
-            u.wx + v.wx,
-            u.xz + v.xz,
-            u.yw + v.yw);
-    }
-
-    Vec mul(float c, Vec v) {
-        return Vec(c * v.x,
-            c * v.y,
-            c * v.z,
-            c * v.w);
-    }
-
-    Bivec mul(float c, Bivec v) {
-        Bivec(c * v.xy,
-            c * v.yz,
-            c * v.zw,
-            c * v.wx,
-            c * v.xz,
-            c * v.yw);
-    }
-
-    Vec mul(Mat m, Vec v) {
-        return add(add(mul(v.x, m.x), mul(v.y, m.y)), add(mul(v.z, m.z), mul(v.w, m.w)));
-    }
-
-    Mat mul(Mat m, Mat n) {
-        return Mat(
-            mul(m, n.x),
-            mul(m, n.y),
-            mul(m, n.z),
-            mul(m, n.w)
-        );
-    }
-
-    Mat mul(float c, Mat m) {
-        return Mat(
-            mul(c, m.x),
-            mul(c, m.y),
-            mul(c, m.z),
-            mul(c, m.w)
-        );
-    }
-
-    Vec normalize(Vec v) {
-        return mul(1 / length(v), v);
-    }
-
-    Bivec normalize(Bivec v) {
-        return mul(1 / length(v), v);
-    }
-
-    Vec tform(float c, Vec v) {
-
-    }
-
-    Vec tform(Vec u, Vec v) {
-
-    }
-
-    Vec tform(Bivec u, Vec v) {
-
-    }
-
-    Mat matrix(float c) {
-        auto m = unit::identity();
-        m.x = tform(c, m.x);
-        m.y = tform(c, m.y);
-        m.z = tform(c, m.z);
-        m.w = tform(c, m.w);
-        return m;
-    }
-
-    Mat matrix(Vec v) {
-        auto m = unit::identity();
-        m.x = tform(v, m.x);
-        m.y = tform(v, m.y);
-        m.z = tform(v, m.z);
-        m.w = tform(v, m.w);
-        return m;
-    }
-
-    Mat matrix(Bivec v) {
-        auto m = unit::identity();
-        m.x = tform(v, m.x);
-        m.y = tform(v, m.y);
-        m.z = tform(v, m.z);
-        m.w = tform(v, m.w);
-        return m;
-    }
-}

main/shaders/main.frag

@@ -1,20 +0,0 @@
-#version 440
-
-#define PI 3.14159
-
-uniform float bright;
-
-in vec4 pos;
-in vec4 pos4;
-in vec4 pos3;
-
-out vec4 color;
-
-void main() {
-    float light = -pos3.y / 10 + .5;
-
-    vec3 col = vec3(bright);
-
-    //color = vec4(col * light, 1);
-    color = vec4(col, 1);
-}

main/shaders/main.geom

@@ -1,47 +0,0 @@
-#version 440
-
-#define CIRCLE_RES 128
-#define PI 3.14159
-
-layout(points) in;
-layout(line_strip, max_vertices=CIRCLE_RES) out;
-
-layout(binding=1) uniform Matrices {
-    mat4 proj;
-    mat4 model;
-};
-
-in float xi_[];
-in float eta_[];
-
-out vec4 pos;
-out vec4 pos4;
-out vec4 pos3;
-
-void main(){
-    for(int k = 0; k <= CIRCLE_RES; k++) {
-        vec2 xi = vec2(xi_[0], 4 * PI * k / (CIRCLE_RES - 1));
-        float eta = eta_[0];
-
-        // todo parameterize the projected circle so there aren't jagged edges.
-        // should be smooth at <=32 segments, not just 128+
-        float x = cos((xi.y + xi.x) / 2) * sin(eta);
-        float y = sin((xi.y + xi.x) / 2) * sin(eta);
-        float z = cos((xi.y - xi.x) / 2) * cos(eta);
-        float w = sin((xi.y - xi.x) / 2) * cos(eta);
-
-        pos = vec4(x, y, z, w);
-        pos4 = model * pos;
-        pos3 = vec4(pos4.xyz / (1 - pos4.w), 1);
-        gl_Position = proj * pos3;
-
-        if (length(gl_Position) > 4) {
-            EndPrimitive();
-            continue;
-        }
-
-        EmitVertex();
-    }
-
-    EndPrimitive();
-}

main/shaders/main.vert

@@ -1,12 +0,0 @@
-#version 440
-
-layout(location=0) in float xi;
-layout(location=1) in float eta;
-
-out float xi_;
-out float eta_;
-
-void main(){
-    xi_ = xi;
-    eta_ = eta;
-}

main/src/main.cpp

@@ -1,397 +0,0 @@
-#include <glad/glad.h>
-#include <GLFW/glfw3.h>
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <vector>
-#include <cmath>
-
-#include "glga.hpp"
-
-#define PI 3.14159f
-
-namespace util {
-    template<typename T>
-    void bufferData(GLenum target, std::vector<T> data, GLenum usage) {
-        glBufferData(target, data.size() * sizeof(T), &data.front(), usage);
-    }
-
-    template<typename T>
-    void bufferData(GLenum target, T &data, GLenum usage) {
-        glBufferData(target, sizeof(T), &data, usage);
-    }
-
-    std::string readFile(const std::string &path) {
-        std::ifstream file(path);
-        if (!file) return std::string();
-
-        file.ignore(std::numeric_limits<std::streamsize>::max());
-        auto size = file.gcount();
-
-        if (size > 0x10000) return std::string();
-
-        file.clear();
-        file.seekg(0, std::ios_base::beg);
-
-        std::stringstream sstr;
-        sstr << file.rdbuf();
-        file.close();
-
-        return sstr.str();
-    }
-
-    void shaderFiles(GLuint shader, std::vector<std::string> &paths) {
-        std::vector<std::string> strs;
-        std::vector<const char *> c_strs;
-
-        for (const auto &path : paths) strs.push_back(readFile(path));
-        for (const auto &str:strs) c_strs.push_back(str.c_str());
-
-        glShaderSource(shader, (GLsizei) c_strs.size(), &c_strs.front(), nullptr);
-    }
-
-    std::string shaderInfoLog(GLuint shader) {
-        GLint log_len;
-        glGetShaderiv(shader, GL_INFO_LOG_LENGTH, &log_len);
-        char log[log_len];
-        glGetShaderInfoLog(shader, log_len, nullptr, log);
-        return std::string(log);
-    }
-
-    std::string programInfoLog(GLuint program) {
-        GLint log_len;
-        glGetProgramiv(program, GL_INFO_LOG_LENGTH, &log_len);
-        char log[log_len];
-        glGetProgramInfoLog(program, log_len, nullptr, log);
-        return std::string(log);
-    }
-
-    GLuint buildShader(GLenum kind, const std::string &name, std::vector<std::string> paths) {
-        GLuint shader = glCreateShader(kind);
-        shaderFiles(shader, paths);
-
-        glCompileShader(shader);
-
-        GLint comp;
-        glGetShaderiv(shader, GL_COMPILE_STATUS, &comp);
-        if (!comp) {
-            std::string log = shaderInfoLog(shader);
-            fprintf(stderr, "SHADER ERROR (%s):\n%s", name.c_str(), log.c_str());
-            glDeleteShader(shader);
-            return 0;
-        }
-
-        return shader;
-    }
-
-    GLuint buildProgram(bool separable, const std::string &name, std::vector<GLuint> shaders) {
-        GLuint program = glCreateProgram();
-
-        if (separable)
-            glProgramParameteri(program, GL_PROGRAM_SEPARABLE, GL_TRUE);
-
-        for (GLuint shader : shaders)
-            glAttachShader(program, shader);
-
-        glLinkProgram(program);
-
-        GLint link;
-        glGetProgramiv(program, GL_LINK_STATUS, &link);
-        if (!link) {
-            std::string log = programInfoLog(program);
-            fprintf(stderr, "PROGRAM ERROR (%s):\n%s", name.c_str(), log.c_str());
-            glDeleteProgram(program);
-            return 0;
-        }
-
-        return program;
-    }
-}
-
-struct HopfCircle {
-    float xi; // longitude
-    float eta; // latitude
-
-    HopfCircle(float xi, float eta) : xi(xi), eta(eta) {}
-
-    HopfCircle() : HopfCircle(0, 0) {}
-};
-
-struct State {
-    GLuint vao{};
-    GLuint vbo{}, ubo{};
-    GLuint unif_bright{};
-    GLuint ubo_bp = 1;
-
-    GLuint prog{};
-
-    std::vector<HopfCircle> circles{};
-    ga::Mat rot = ga::unit::identity();
-    float view = 4.f;
-
-    float t = 0;
-    float dt = 0;
-
-    void regen() {
-        circles.clear();
-
-        const int N = 32;
-        const int M = 4;
-        const float PAD = 0.0125;
-
-        for (int x = 0; x <= N; ++x) {
-            for (int e = 1; e < M; ++e) {
-                float xi = 2 * PI * x / N;
-                float eta = (PI / 2 - 2 * PAD) * e / M + PAD;
-                circles.emplace_back(xi, eta);
-            }
-        }
-
-        glBindBuffer(GL_ARRAY_BUFFER, vbo);
-        util::bufferData(GL_ARRAY_BUFFER, circles, GL_STREAM_DRAW);
-        glBindBuffer(GL_ARRAY_BUFFER, 0);
-    }
-
-    explicit State(GLFWwindow *window) {
-        printf("vendor: %s\nrenderer: %s\n", glGetString(GL_VENDOR), glGetString(GL_RENDERER));
-
-        float lw[2];
-        glGetFloatv(GL_LINE_WIDTH_RANGE, lw);
-        printf("line width range: %.2f to %.2f\n", lw[0], lw[1]);
-
-        glGenBuffers(1, &vbo);
-        glGenBuffers(1, &ubo);
-        glGenVertexArrays(1, &vao);
-
-        regen();
-
-        glBindBufferBase(GL_UNIFORM_BUFFER, ubo_bp, ubo);
-
-        GLuint vs = util::buildShader(GL_VERTEX_SHADER, "vs", {"shaders/main.vert"});
-        GLuint gs = util::buildShader(GL_GEOMETRY_SHADER, "gs", {"shaders/main.geom"});
-        GLuint fs = util::buildShader(GL_FRAGMENT_SHADER, "fs", {"shaders/main.frag"});
-        prog = util::buildProgram(false, "prog", {vs, gs, fs});
-
-        unif_bright = (GLuint) glGetUniformLocation(prog, "bright");
-
-        glBindVertexArray(vao);
-
-        glBindBuffer(GL_ARRAY_BUFFER, vbo);
-        auto xi = glGetAttribLocation(prog, "xi");
-        if (xi >= 0) {
-            glEnableVertexAttribArray((unsigned) xi);
-            glVertexAttribPointer((unsigned) xi, 1, GL_FLOAT, GL_FALSE,
-                sizeof(HopfCircle),
-                (void *) offsetof(HopfCircle, xi));
-        }
-
-        auto eta = glGetAttribLocation(prog, "eta");
-        if (eta >= 0) {
-            glEnableVertexAttribArray((unsigned) eta);
-            glVertexAttribPointer((unsigned) eta, 1, GL_FLOAT, GL_FALSE,
-                sizeof(HopfCircle),
-                (void *) offsetof(HopfCircle, eta));
-        }
-        glBindBuffer(GL_ARRAY_BUFFER, 0);
-
-        glBindVertexArray(0);
-    }
-
-    void update(GLFWwindow *window, float dt, int frame) {
-        int w, h;
-        glfwGetFramebufferSize(window, &w, &h);
-        auto ar = (float) w / h;
-
-        t += dt;
-
-        float c = std::cos(t / 8);
-        float s = std::sin(t / 8);
-
-        float c_ = std::cos(3.f / 8);
-        float s_ = std::sin(3.f / 8);
-
-        glBindBuffer(GL_UNIFORM_BUFFER, ubo);
-        util::bufferData<ga::Mat>(GL_UNIFORM_BUFFER, {
-            ga::Mat(
-                ga::Vec(1.f / view / ar, 0, 0, 0),
-                ga::Vec(0, 0, 1.f / 10, 0),
-                ga::Vec(0, 1.f / view, 0, 0),
-                ga::Vec(0, 0, 0, 1.f)),
-            rot
-        }, GL_STREAM_DRAW);
-        glBindBuffer(GL_UNIFORM_BUFFER, 0);
-    }
-
-    void render(GLFWwindow *window, float dt, int frame) {
-        int w, h;
-        glfwGetFramebufferSize(window, &w, &h);
-        glViewport(0, 0, w, h);
-
-        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-
-        glDepthFunc(GL_LEQUAL);
-
-        glUseProgram(prog);
-        glBindVertexArray(vao);
-        glEnable(GL_DEPTH_TEST);
-
-        glLineWidth(10);
-        glUniform1f(unif_bright, 0);
-        glDrawArrays(GL_POINTS, 0, (unsigned) circles.size());
-
-        glLineWidth(1);
-        glUniform1f(unif_bright, 1);
-        glDrawArrays(GL_POINTS, 0, (unsigned) circles.size());
-
-        glBindVertexArray(0);
-        glUseProgram(0);
-
-        glFinish();
-        glfwSwapBuffers(window);
-    }
-
-    void on_scroll(GLFWwindow *window, double xoffset, double yoffset) {
-        printf("scroll %.2f %.2f\n", xoffset, yoffset);
-
-        int xw = glfwGetKey(window, GLFW_KEY_Q);
-        int yw = glfwGetKey(window, GLFW_KEY_W);
-        int zw = glfwGetKey(window, GLFW_KEY_E);
-        int xy = glfwGetKey(window, GLFW_KEY_A);
-        int yz = glfwGetKey(window, GLFW_KEY_S);
-        int zx = glfwGetKey(window, GLFW_KEY_D);
-
-        int zoom = glfwGetKey(window, GLFW_KEY_LEFT_CONTROL);
-
-        float t = (float) yoffset * PI / 48;
-        float c = std::cos(t);
-        float s = std::sin(t);
-
-        if (xw) {
-            rot = ga::mul(rot, ga::Mat(
-                ga::Vec(c, 0, 0, s),
-                ga::Vec(0, 1, 0, 0),
-                ga::Vec(0, 0, 1, 0),
-                ga::Vec(-s, 0, 0, c)
-            ));
-        }
-
-        if (yw) {
-            rot = ga::mul(rot, ga::Mat(
-                ga::Vec(1, 0, 0, 0),
-                ga::Vec(0, c, 0, s),
-                ga::Vec(0, 0, 1, 0),
-                ga::Vec(0, -s, 0, c)
-            ));
-        }
-
-        if (zw) {
-            rot = ga::mul(rot, ga::Mat(
-                ga::Vec(1, 0, 0, 0),
-                ga::Vec(0, 1, 0, 0),
-                ga::Vec(0, 0, c, s),
-                ga::Vec(0, 0, -s, c)
-            ));
-        }
-
-        if (xy) {
-            rot = ga::mul(rot, ga::Mat(
-                ga::Vec(c, s, 0, 0),
-                ga::Vec(-s, c, 0, 0),
-                ga::Vec(0, 0, 1, 0),
-                ga::Vec(0, 0, 0, 1)
-            ));
-        }
-
-        if (yz) {
-            rot = ga::mul(rot, ga::Mat(
-                ga::Vec(1, 0, 0, 0),
-                ga::Vec(0, c, s, 0),
-                ga::Vec(0, -s, c, 0),
-                ga::Vec(0, 0, 0, 1)
-            ));
-        }
-
-        if (zx) {
-            rot = ga::mul(rot, ga::Mat(
-                ga::Vec(c, 0, -s, 0),
-                ga::Vec(0, 1, 0, 0),
-                ga::Vec(s, 0, c, 0),
-                ga::Vec(0, 0, 0, 1)
-            ));
-        }
-
-        if (zoom) {
-            view *= 1 - t / 5;
-        }
-
-        std::ofstream matfile;
-        matfile.open("matrix.txt");
-        matfile << "vec((\n";
-        matfile << "vec((" << rot.x.x << ", " << rot.x.y << ", " << rot.x.z << ", " << rot.x.w << ")),\n";
-        matfile << "vec((" << rot.y.x << ", " << rot.y.y << ", " << rot.y.z << ", " << rot.y.w << ")),\n";
-        matfile << "vec((" << rot.z.x << ", " << rot.z.y << ", " << rot.z.z << ", " << rot.z.w << ")),\n";
-        matfile << "vec((" << rot.w.x << ", " << rot.w.y << ", " << rot.w.z << ", " << rot.w.w << ")),\n";
-        matfile << ")),";
-        matfile.close();
-    }
-
-    void deinit(GLFWwindow *window) {
-    }
-};
-
-void scroll_callback(GLFWwindow *window, double xoffset, double yoffset) {
-    auto *state = (State *) glfwGetWindowUserPointer(window);
-    state->on_scroll(window, xoffset, yoffset);
-}
-
-void run() {
-    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4);
-    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 4);
-    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
-
-    auto window = glfwCreateWindow(1280, 720, "Hopf Fibration", nullptr, nullptr);
-    if (!window) {
-        glfwTerminate();
-        exit(EXIT_FAILURE);
-    }
-
-    glfwMakeContextCurrent(window);
-    gladLoadGLLoader((GLADloadproc) glfwGetProcAddress);
-    glfwSwapInterval(0);
-
-    auto state = State(window);
-    glfwSetWindowUserPointer(window, &state);
-
-    glfwSetScrollCallback(window, scroll_callback);
-
-    double time = glfwGetTime();
-    int frame = 0;
-    while (!glfwWindowShouldClose(window)) {
-        double time_ = glfwGetTime();
-        auto dt = (float) (time_ - time);
-
-        state.update(window, dt, frame);
-        state.render(window, dt, frame);
-
-        glfwPollEvents();
-
-        time = time_;
-        frame += 1;
-    }
-
-    state.deinit(window);
-
-    glfwDestroyWindow(window);
-}
-
-int main() {
-    if (!glfwInit()) {
-        return EXIT_FAILURE;
-    }
-
-    run();
-
-    glfwTerminate();
-    return EXIT_SUCCESS;
-}

res/shaders/4d.vert.glsl

@@ -0,0 +1,15 @@
+#version 440
+
+layout(location=1) uniform float time;
+layout(location=2) uniform mat4 proj;
+layout(location=3) uniform mat4 rot;
+layout(location=4) uniform mat4 view;
+
+layout(location=0) in vec4 pos;
+
+layout(location=0) out vec4 opos;
+
+void main() {
+    opos = rot * pos;
+    gl_Position = proj * view * vec4(opos.xyz / (1 - opos.w), 1.0);
+}

res/shaders/main-discard.frag.glsl

@@ -0,0 +1,17 @@
+#version 440
+
+layout(location=0) uniform vec4 ucol;
+
+layout(location=0) in vec4 pos;
+
+layout(location=0) out vec4 col;
+
+void main() {
+    vec3 p = pos.xyz / (1 - pos.w);
+    if (dot(p, p) > 3) discard;
+    
+    float d = 1.0 - gl_FragCoord.z;
+    d = (d - 0.5) / 0.7 + 0.5;
+    col = ucol;
+    col.xyz *= d;
+}

res/shaders/main.frag.glsl

@@ -0,0 +1,12 @@
+#version 440
+
+layout(location=0) uniform vec4 ucol;
+
+layout(location=0) out vec4 col;
+
+void main() {
+    float d = 1.0 - gl_FragCoord.z;
+    d = (d - 0.5) / 0.7 + 0.5;
+    col = ucol;
+    col.xyz *= d;
+}

res/shaders/main.vert.glsl

@@ -0,0 +1,13 @@
+#version 440
+
+layout(location=1) uniform float time;
+layout(location=2) uniform mat4 proj;
+layout(location=3) uniform mat4 rot;
+layout(location=0) in vec3 pos;
+
+void main() {
+    mat3 rot3 = mat3(rot);
+    vec3 pos3 = rot3 * pos;
+
+    gl_Position = proj * vec4(pos3, 1.0);
+}

src/gl/buffer.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <utility>
+
+template<typename T_>
+class Buffer {
+public:
+    using Type = T_;
+
+private:
+    GLuint id = 0;
+
+public:
+    Buffer() {
+        glCreateBuffers(1, &id);
+    }
+
+    Buffer(Buffer&& o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    Buffer(const Buffer&) = delete;  // this is doable, but would be slow.
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+
+    template<typename RandomIt>
+    GLuint upload(RandomIt first, RandomIt last, GLenum mode = GL_STATIC_DRAW) {
+        size_t count = last - first;
+
+        // todo StaticBuffer that uses BufferStorage
+        glNamedBufferData(id, sizeof(Type) * count, nullptr, mode);
+
+        Type* out = (Type*) glMapNamedBuffer(id, GL_WRITE_ONLY);
+        std::copy(first, last, out);
+        glUnmapNamedBuffer(id);
+
+        return count;
+    }
+
+    template<typename T>
+    GLuint upload(const T& data, GLenum mode = GL_STATIC_DRAW) {
+        return upload(data.begin(), data.end(), mode);
+    }
+
+    ~Buffer() {
+        // delete silently ignores 0.
+        glDeleteBuffers(1, &id);
+    }
+};

src/gl/debug.hpp

@@ -0,0 +1,80 @@
+#pragma once
+
+#ifndef NDEBUG
+
+void GLAPIENTRY log_gl_debug_callback(
+    GLenum source,
+    GLenum type,
+    GLuint id,
+    GLenum severity,
+    GLsizei length,
+    const GLchar *message,
+    const void *userParam
+) {
+    std::string s_source;
+    switch(type){
+        case GL_DEBUG_SOURCE_API:
+            s_source = "API:";
+        case GL_DEBUG_SOURCE_WINDOW_SYSTEM:
+            s_source = "WINDOW:";
+        case GL_DEBUG_SOURCE_SHADER_COMPILER:
+            s_source = "SHADER:";
+        case GL_DEBUG_SOURCE_THIRD_PARTY:
+            s_source = "3P:";
+        case GL_DEBUG_SOURCE_APPLICATION:
+            s_source = "APP:";
+        default:
+            s_source = "";
+    }
+
+    std::string s_type;
+    switch (type) {
+    case GL_DEBUG_TYPE_ERROR:
+        s_type = "ERROR:";
+        break;
+    case GL_DEBUG_TYPE_DEPRECATED_BEHAVIOR:
+        s_type = "DEPRECATED:";
+        break;
+    case GL_DEBUG_TYPE_UNDEFINED_BEHAVIOR:
+        s_type = "UNDEFINED:";
+        break;
+    case GL_DEBUG_TYPE_PORTABILITY:
+        s_type = "PORTABILITY:";
+        break;
+    case GL_DEBUG_TYPE_PERFORMANCE:
+        s_type = "PERFORMANCE:";
+        break;
+    case GL_DEBUG_TYPE_MARKER:
+        s_type = "MARKER:";
+        break;
+    case GL_DEBUG_TYPE_PUSH_GROUP:
+        s_type = "PUSH_GROUP:";
+        break;
+    case GL_DEBUG_TYPE_POP_GROUP:
+        s_type = "POP_GROUP:";
+        break;
+    default:
+        s_type = "";
+        break;
+    }
+
+    std::string s_severity;
+    switch (severity) {
+    case GL_DEBUG_SEVERITY_HIGH:
+        s_severity = "HIGH:";
+        break;
+    case GL_DEBUG_SEVERITY_MEDIUM:
+        s_severity = "MED:";
+        break;
+    case GL_DEBUG_SEVERITY_LOW:
+        s_severity = "LOW:";
+        break;
+    default:
+        s_severity = "INFO:";
+        break;
+    }
+
+    std::cerr << "GL:" << s_source << s_type << s_severity << " " << message << std::endl;
+}
+
+#endif

src/gl/shader.hpp

@@ -0,0 +1,79 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <string>
+#include <memory>
+#include <vector>
+
+template<GLenum mode>
+class Shader {
+public:
+private:
+    GLuint id = 0;
+
+public:
+    explicit Shader(const std::string &src) {
+        id = glCreateShader(mode);
+
+        const char *str = src.c_str();
+        glShaderSource(id, 1, &str, nullptr);
+        glCompileShader(id);
+
+        // todo throw if compile failed
+    }
+
+    explicit Shader(std::ifstream source)
+        : Shader(std::string(
+        std::istreambuf_iterator<char>(source),
+        std::istreambuf_iterator<char>()
+    )) {}
+
+    Shader(const Shader &) = delete;
+
+    Shader(Shader &&o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    ~Shader() {
+        glDeleteShader(id);
+    }
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+};
+
+using VertexShader = Shader<GL_VERTEX_SHADER>;
+using FragmentShader = Shader<GL_FRAGMENT_SHADER>;
+
+class Program {
+private:
+    GLuint id = 0;
+
+public:
+    template<GLenum ...mode>
+    explicit Program(const Shader<mode> &...shader) {
+        id = glCreateProgram();
+
+        (glAttachShader(id, shader), ...);
+
+        glLinkProgram(id);
+
+        // todo throw if link failed
+    }
+
+    Program(const Program &) = delete;
+
+    Program(Program &&o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    ~Program() {
+        glDeleteProgram(id);
+    }
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+};

src/gl/types.hpp

@@ -0,0 +1,118 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <Eigen/Eigen>
+
+template<typename T_>
+struct Buffer;
+
+template<class ...Fmt_>
+struct VertexArray;
+
+
+template<typename ctype, GLenum type, GLint size, GLuint offset = 0, GLboolean norm = GL_FALSE>
+struct Format {
+    using CType = ctype;
+    static constexpr GLenum Type = type;
+    static constexpr GLint Size = size;
+    static constexpr GLuint Offset = offset;
+    static constexpr GLboolean Norm = norm;
+
+    template<typename ctype_, GLint size_>
+    using As = Format<ctype_, type, size_, offset, norm>;
+
+    static inline constexpr void apply(GLuint vao, GLuint idx) {
+        glEnableVertexArrayAttrib(vao, idx);
+        glVertexArrayAttribFormat(vao, idx, size, type, norm, offset);
+    }
+};
+
+template<typename ctype, GLenum type, GLint size, GLuint offset = 0>
+struct IFormat {
+    using CType = ctype;
+    static constexpr GLenum Type = type;
+    static constexpr GLint Size = size;
+    static constexpr GLuint Offset = offset;
+
+    template<typename ctype_, GLint size_>
+    using As = IFormat<ctype_, type, size_, offset>;
+
+    static inline constexpr void apply(GLuint vao, GLuint idx) {
+        glEnableVertexArrayAttrib(vao, idx);
+        glVertexArrayAttribIFormat(vao, idx, size, type, offset);
+    }
+};
+
+template<typename ctype, GLenum type, GLint size, GLuint offset = 0>
+struct LFormat {
+    using CType = ctype;
+    static constexpr GLenum Type = type;
+    static constexpr GLint Size = size;
+    static constexpr GLuint Offset = offset;
+
+    template<typename ctype_, GLint size_>
+    using As = LFormat<ctype_, type, size_, offset>;
+
+    static inline constexpr void apply(GLuint vao, GLuint idx) {
+        glEnableVertexArrayAttrib(vao, idx);
+        glVertexArrayAttribLFormat(vao, idx, size, type, offset);
+    }
+};
+
+template<typename Fmt_>
+struct AutoFormat {
+    using Fmt = Fmt_;
+};
+
+template<>
+struct AutoFormat<GLdouble> {
+    using Fmt = LFormat<GLdouble, GL_DOUBLE, 1>;
+};
+
+template<>
+struct AutoFormat<GLfloat> {
+    using Fmt = Format<GLfloat, GL_FLOAT, 1>;
+};
+
+template<>
+struct AutoFormat<GLint> {
+    using Fmt = IFormat<GLint, GL_INT, 1>;
+};
+
+template<>
+struct AutoFormat<GLuint> {
+    using Fmt = IFormat<GLuint, GL_UNSIGNED_INT, 1>;
+};
+
+template<typename Scalar, int Dim>
+struct AutoFormat<Eigen::Vector<Scalar, Dim>> {
+    static_assert((Dim >= 1) && (Dim <= 4), "Dim must be in range 1-4");
+
+    using Fmt = typename AutoFormat<Scalar>::Fmt::template As<Eigen::Vector<Scalar, Dim>, Dim>;
+};
+
+template<typename Fmt_>
+struct Binder {
+    using CType = typename AutoFormat<Fmt_>::Fmt::CType;
+
+    const GLuint buf;
+    const GLuint offset;
+    const GLuint stride;
+
+    Binder(const Buffer<CType> &buf) // NOLINT(google-explicit-constructor)
+        : buf((GLuint) buf), offset(0), stride(sizeof(CType)) {}
+
+    template<typename CType_>
+    Binder(const Buffer<CType_> &buf, GLuint offset)
+        : buf((GLuint) buf), offset(offset), stride(sizeof(CType_)) {}
+
+    Binder(const Binder<Fmt_> &o)
+        : buf(o.buf), offset(o.offset), stride(o.stride) {}
+
+    inline void bind(GLuint vao, GLuint idx) const {
+        glVertexArrayVertexBuffer(vao, idx, buf, offset, stride);
+    }
+};
+
+#define ATTR(buf, field) Binder<decltype(decltype(buf)::Type::field)>(buf, offsetof(decltype(buf)::Type, field))

src/gl/vertexarray.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <utility>
+
+#include "buffer.hpp"
+#include "types.hpp"
+
+template<typename ...Fmt_>
+class VertexArray {
+public:
+    template<size_t idx>
+    using Fmt = std::tuple_element_t<idx, std::tuple<Fmt_...>>;
+
+private:
+    GLuint id = 0;
+
+    template<size_t ...idx>
+    inline void formatall(
+        std::integer_sequence<size_t, idx...>
+    ) {
+        (AutoFormat<Fmt_>::Fmt::apply(id, idx), ...);
+    }
+
+    template<size_t ...idx>
+    inline void bindall(
+        const Binder<Fmt_> &...buf,
+        std::integer_sequence<size_t, idx...>
+    ) {
+        (buf.bind(id, idx), ...);
+    }
+
+public:
+    explicit VertexArray() {
+        glCreateVertexArrays(1, &id);
+        formatall(std::make_index_sequence<sizeof...(Fmt_)>());
+    }
+
+    explicit VertexArray(const Binder<Fmt_> &...buf)
+        : VertexArray() {
+        bind(buf...);
+    }
+
+    VertexArray(VertexArray &&o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    VertexArray(const VertexArray &) = delete;  // this is doable, but would be slow.
+
+    ~VertexArray() {
+        glDeleteVertexArrays(1, &id);
+    }
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+
+    void bind(const Binder<Fmt_> &...buf) {
+        bindall(buf..., std::make_index_sequence<sizeof...(Fmt_)>());
+    }
+
+    template<size_t idx>
+    void bind(const Binder<Fmt<idx>> &buf) {
+        buf.bind(id, idx);
+    }
+};

src/main.cpp

@@ -0,0 +1,325 @@
+#include <glad/glad.h>
+#include <GLFW/glfw3.h>
+#include <imgui.h>
+#include <backends/imgui_impl_glfw.h>
+#include <backends/imgui_impl_opengl3.h>
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+#include "gl/debug.hpp"
+#include "gl/buffer.hpp"
+#include "gl/shader.hpp"
+#include "gl/vertexarray.hpp"
+
+#include <ml/meshlib.hpp>
+#include <ml/meshlib_json.hpp>
+
+struct State {
+    using Transform4 = Eigen::Transform<float, 4, Eigen::Projective>;
+    using Transform3 = Eigen::Transform<float, 3, Eigen::Projective>;
+
+    Eigen::Vector4f bg{0.07f, 0.09f, 0.10f, 1.00f};
+    Eigen::Vector4f fg{0.71f, 0.53f, 0.94f, 1.00f};
+    Eigen::Vector4f wf{0.95f, 0.95f, 0.95f, 1.00f};
+
+    Transform4 rot = Transform4::Identity();
+    Transform3 view = Transform3::Identity();
+};
+
+Eigen::Matrix4f rotor(int u, int v, float rad) {
+    Eigen::Matrix4f res = Eigen::Matrix4f::Identity();
+    res(u, u) = res(v, v) = cosf(rad);
+    res(u, v) = res(v, u) = sinf(rad);
+    res(u, v) *= -1;
+    return res;
+}
+
+template<typename T_>
+T_ mix(const T_ &a, const T_ &b, const typename T_::Scalar &x) {
+    return a * (1 - x) + b * x;
+}
+
+void show_overlay(State &state) {
+    static std::string gl_vendor = (const char *) glGetString(GL_VENDOR);
+    static std::string gl_renderer = (const char *) glGetString(GL_RENDERER);
+    static std::string gl_version = (const char *) glGetString(GL_VERSION);
+    static std::string glsl_version = (const char *) glGetString(GL_SHADING_LANGUAGE_VERSION);
+
+    ImGuiWindowFlags window_flags =
+        ImGuiWindowFlags_AlwaysAutoResize |
+            ImGuiWindowFlags_NoSavedSettings |
+            ImGuiWindowFlags_NoFocusOnAppearing |
+            ImGuiWindowFlags_NoNav |
+            ImGuiWindowFlags_NoBringToFrontOnFocus |
+            ImGuiWindowFlags_NoMove;
+
+    ImGuiStyle &style = ImGui::GetStyle();
+    const auto PAD = style.DisplaySafeAreaPadding;
+    auto window_pos = PAD;
+
+    ImGui::SetNextWindowPos(window_pos, ImGuiCond_Always);
+    ImGui::SetNextWindowBgAlpha(0.35f * style.Alpha);
+    ImGui::SetNextWindowCollapsed(true, ImGuiCond_Appearing);
+    ImGui::Begin("Graphics Information", nullptr, window_flags);
+    ImGui::Text("GL Vendor    | %s", gl_vendor.c_str());
+    ImGui::Text("GL Renderer  | %s", gl_renderer.c_str());
+    ImGui::Text("GL Version   | %s", gl_version.c_str());
+    ImGui::Text("GLSL Version | %s", glsl_version.c_str());
+
+    auto v2 = ImGui::GetWindowSize();
+    window_pos.y += v2.y + PAD.y;
+    ImGui::End();
+
+    ImGui::SetNextWindowPos(window_pos, ImGuiCond_Always);
+    ImGui::SetNextWindowBgAlpha(0.35f * style.Alpha);
+    ImGui::SetNextWindowCollapsed(true, ImGuiCond_Appearing);
+    ImGui::Begin("Controls", nullptr, window_flags);
+
+    ImGuiIO &io = ImGui::GetIO();
+    ImGui::Text("FPS          | %.2f", io.Framerate);
+
+    ImGui::Separator();
+
+    ImGui::ColorEdit3("Background", state.bg.data(), ImGuiColorEditFlags_Float);
+    ImGui::ColorEdit3("Foreground", state.fg.data(), ImGuiColorEditFlags_Float);
+    ImGui::ColorEdit3("Wireframe", state.wf.data(), ImGuiColorEditFlags_Float);
+
+    if (io.MouseDown[0] && !io.WantCaptureMouse) {
+        Eigen::Matrix4f rot = Eigen::Matrix4f::Identity();
+        Eigen::Vector2f del{io.MouseDelta.x, io.MouseDelta.y};
+        del /= 200.0f;
+
+        if (io.KeyShift) {
+            del /= 5.0f;
+        }
+
+        if (io.KeyCtrl) {
+            Eigen::Matrix4f rx = rotor(0, 3, -del.x());
+            Eigen::Matrix4f ry = rotor(1, 3, del.y());
+            rot = rx * ry;
+        } else {
+            Eigen::Matrix4f rx = rotor(0, 2, -del.x());
+            Eigen::Matrix4f ry = rotor(1, 2, del.y());
+            rot = rx * ry;
+        }
+
+        state.rot.linear() = rot * state.rot.linear();
+    }
+
+    if ((io.MouseWheel != 0) && !io.WantCaptureMouse) {
+        float scale = 1.0f + io.MouseWheel * 0.05f;
+        state.view.linear() *= scale;
+    }
+
+    ImGui::End();
+}
+
+void set_style() {
+    ImGui::StyleColorsDark();
+
+    ImGuiStyle &style = ImGui::GetStyle();
+    style.WindowRounding = 4;
+    style.FrameRounding = 2;
+    style.DisplaySafeAreaPadding.x = 10;
+    style.DisplaySafeAreaPadding.y = 10;
+}
+
+
+Eigen::Vector4f hopf_map(float e1, float e2, float n) {
+    Eigen::Vector4f res;
+    res <<
+        cosf((e2 + e1) / 2) * sinf(n),
+        sinf((e2 + e1) / 2) * sinf(n),
+        cosf((e2 - e1) / 2) * cosf(n),
+        sinf((e2 - e1) / 2) * cosf(n);
+    return res;
+}
+
+auto make_hopf(size_t latitudes, size_t longitudes, size_t link_res) {
+    Eigen::Matrix4Xf points(4, latitudes * longitudes * link_res);
+    ml::Matrix1Xui lines(1, latitudes * longitudes * link_res);
+
+    Eigen::Index idx = 0;
+
+    for (int i = 0; i < latitudes; ++i) {
+        float n = (float) i / (float) (latitudes - 1) * M_PIf32 / 2.0f;
+
+        for (int j = 0; j < longitudes; ++j) {
+            float e1 = (float) j / (float) longitudes * M_PIf32 * 2.0f;
+
+            for (int k = 0; k < link_res; ++k) {
+                float e2 = (float) k / (float) link_res * M_PIf32 * 4.0f;
+
+                lines.col(idx) << idx;
+                points.col(idx) = hopf_map(e1, e2, n).normalized();
+                idx++;
+            }
+        }
+    }
+
+    return ml::Mesh(points, lines);
+}
+
+auto make_link(float n, float e1, size_t link_res) {
+    Eigen::Matrix4Xf points(4, link_res);
+    ml::Matrix1Xui lines(1, link_res);
+
+    Eigen::Index idx = 0;
+
+    for (int k = 0; k < link_res; ++k) {
+        float e2 = (float) k / (float) link_res * M_PIf32 * 4.0f;
+
+        lines.col(idx) << idx;
+        points.col(idx) = hopf_map(e1, e2, n).normalized();
+        idx++;
+    }
+
+    return ml::Mesh(points, lines);
+}
+
+int run(GLFWwindow *window, ImGuiContext *context) {
+    State state;
+
+    state.fg *= 0.5;
+
+    Buffer<GLuint> ind_buf;
+    Buffer<Eigen::Vector4f> vert_buf;
+    Buffer<Eigen::Vector4f> wire_buf;
+
+    VertexArray<Eigen::Vector4f> vao(vert_buf);
+    glVertexArrayElementBuffer(vao, ind_buf);
+
+    VertexArray<Eigen::Vector4f> wire_vao(wire_buf);
+
+    using PointsType = Eigen::Matrix<float, 4, Eigen::Dynamic>;
+    using CellsType = Eigen::Matrix<unsigned, 3, Eigen::Dynamic>;
+    using Mesh = ml::Mesh<PointsType, CellsType>;
+
+    auto mesh = make_hopf(24, 48, 1024);
+//    auto mesh = make_hopf(5, 32, 1024);
+
+//        auto wire = make_link(0.0f, 0.0f, 256);
+    auto wire = make_hopf(4, 48, 1024);
+    wire_buf.upload(wire.points.colwise(), GL_STREAM_DRAW);
+    
+    auto elements = (GLint) ind_buf.upload(mesh.cells.reshaped());
+    vert_buf.upload(mesh.points.colwise());
+
+    VertexShader vs4d(std::ifstream("res/shaders/4d.vert.glsl"));
+    FragmentShader fs(std::ifstream("res/shaders/main.frag.glsl"));
+    FragmentShader fsd(std::ifstream("res/shaders/main-discard.frag.glsl"));
+
+    Program pgm(vs4d, fs);
+    Program pgm_discard(vs4d, fsd);
+
+    glEnable(GL_DEPTH_TEST);
+
+    Eigen::Projective3f proj;
+
+    while (!glfwWindowShouldClose(window)) {
+        glfwPollEvents();
+
+        ImGui_ImplOpenGL3_NewFrame();
+        ImGui_ImplGlfw_NewFrame();
+        ImGui::NewFrame();
+        show_overlay(state);
+        ImGui::Render();
+
+        int display_w, display_h;
+        glfwGetFramebufferSize(window, &display_w, &display_h);
+        glViewport(0, 0, display_w, display_h);
+        glClearColor(state.bg[0], state.bg[1], state.bg[2], state.bg[3]);
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+        auto aspect = (float) display_h / (float) display_w;
+        proj = Eigen::AlignedScaling3f(aspect, 1.0, -0.6);
+
+        Eigen::Matrix4f rot = state.rot.linear();
+        Eigen::Matrix4f view = state.view.matrix();
+
+        glUseProgram(pgm);
+        glBindVertexArray(vao);
+        glUniform4fv(0, 1, state.fg.data());
+        glUniform1f(1, (GLfloat) glfwGetTime());
+        glUniformMatrix4fv(2, 1, false, proj.data());
+        glUniformMatrix4fv(3, 1, false, rot.data());
+        glUniformMatrix4fv(4, 1, false, view.data());
+        glDrawElements(GL_POINTS, elements, GL_UNSIGNED_INT, nullptr);
+        glBindVertexArray(0);
+        glUseProgram(0);
+
+        glClear(GL_DEPTH_BUFFER_BIT);
+
+        glUseProgram(pgm_discard);
+        glBindVertexArray(wire_vao);
+        glUniform4fv(0, 1, state.wf.data());
+        glUniform1f(1, (GLfloat) glfwGetTime());
+        glUniformMatrix4fv(2, 1, false, proj.data());
+        glUniformMatrix4fv(3, 1, false, rot.data());
+        glUniformMatrix4fv(4, 1, false, view.data());
+        glDrawArrays(GL_POINTS, 0, (GLint) wire.points.cols());
+        glBindVertexArray(0);
+        glUseProgram(0);
+
+        ImGui_ImplOpenGL3_RenderDrawData(ImGui::GetDrawData());
+        glfwSwapBuffers(window);
+    }
+
+    return EXIT_SUCCESS;
+}
+
+int main() {
+    if (!glfwInit()) {
+        std::cerr << "GLFW:Failed initialization" << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4);
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 6);
+    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
+    auto *window = glfwCreateWindow(1280, 720, "Hopf Fibration", nullptr, nullptr);
+    if (!window) {
+        std::cerr << "GLFW:Failed to create window" << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    glfwMakeContextCurrent(window);
+    glfwSwapInterval(1);
+    gladLoadGLLoader((GLADloadproc) glfwGetProcAddress);
+
+#ifndef NDEBUG
+    glEnable(GL_DEBUG_OUTPUT);
+    glDebugMessageCallback(log_gl_debug_callback, nullptr);
+    glDebugMessageControl(
+        GL_DONT_CARE, GL_DEBUG_TYPE_OTHER,
+        GL_DEBUG_SEVERITY_NOTIFICATION,
+        0, nullptr, GL_FALSE
+    );
+#endif
+
+    IMGUI_CHECKVERSION();
+    auto *context = ImGui::CreateContext();
+    ImGui_ImplGlfw_InitForOpenGL(window, true);
+    ImGui_ImplOpenGL3_Init("#version 130");
+
+    set_style();
+
+    int exit_code = EXIT_SUCCESS;
+
+    try {
+        exit_code = run(window, context);
+    } catch (const std::exception &e) {
+        std::cerr << e.what() << std::endl;
+        exit_code = EXIT_FAILURE;
+    }
+
+    ImGui_ImplOpenGL3_Shutdown();
+    ImGui_ImplGlfw_Shutdown();
+    ImGui::DestroyContext();
+
+    glfwDestroyWindow(window);
+    glfwTerminate();
+
+    return exit_code;
+}