Multiple passes for different point-cloud resolutions

.gitignore

@@ -1,150 +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

.gitmodules

@@ -1,10 +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
-[submodule "vendor/glm"]
-	path = vendor/glm
-	url = https://github.com/g-truc/glm.git

.idea/encodings.xml

@@ -1,4 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project version="4">
-  <component name="Encoding" addBOMForNewFiles="with NO BOM" />
-</project>

.idea/hopf-fibration.iml

@@ -1,2 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<module classpath="CMake" type="CPP_MODULE" version="4" />

.idea/misc.xml

@@ -1,12 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project version="4">
-  <component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$" />
-  <component name="CidrRootsConfiguration">
-    <sourceRoots>
-      <file path="$PROJECT_DIR$/main" />
-    </sourceRoots>
-    <libraryRoots>
-      <file path="$PROJECT_DIR$/vendor" />
-    </libraryRoots>
-  </component>
-</project>

.idea/modules.xml

@@ -1,8 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project version="4">
-  <component name="ProjectModuleManager">
-    <modules>
-      <module fileurl="file://$PROJECT_DIR$/.idea/hopf-fibration.iml" filepath="$PROJECT_DIR$/.idea/hopf-fibration.iml" />
-    </modules>
-  </component>
-</project>

.idea/vcs.xml

@@ -1,6 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project version="4">
-  <component name="VcsDirectoryMappings">
-    <mapping directory="$PROJECT_DIR$" vcs="Git" />
-  </component>
-</project>

CMakeLists.txt

@@ -1,17 +1,31 @@
-cmake_minimum_required(VERSION 3.10)
-project(hopf-fibration)
+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)
 
-option(GLM_TEST_ENABLE OFF)
-add_subdirectory(vendor/glm)
+include_directories(include)
 
-add_subdirectory(main)
+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)

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,30 +0,0 @@
-project(main)
-
-add_executable(${PROJECT_NAME}
-        src/main.cpp include/util.h)
-
-target_link_libraries(${PROJECT_NAME}
-        glad
-        glm
-        glfw)
-
-target_include_directories(${PROJECT_NAME}
-        PRIVATE
-        include)
-
-set(SHADERS
-        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/util.h

@@ -1,113 +0,0 @@
-//
-// Created by allem on 2/7/2019.
-//
-
-#ifndef HOPF_FIBRATION_UTIL_H
-#define HOPF_FIBRATION_UTIL_H
-
-#include <glad/glad.h>
-
-#include <fstream>
-#include <sstream>
-#include <string>
-#include <vector>
-
-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;
-    }
-}
-
-#endif //HOPF_FIBRATION_UTIL_H

main/shaders/main.frag

@@ -1,12 +0,0 @@
-#version 440
-
-layout(binding=1) uniform Unifs {
-    mat4 uProj;
-    vec4 uColor;
-};
-
-out vec4 color;
-
-void main() {
-    color = uColor;
-}

main/shaders/main.vert

@@ -1,12 +0,0 @@
-#version 440
-
-layout(binding=1) uniform Unifs {
-    mat4 uProj;
-    vec4 uColor;
-};
-
-in vec4 iPos;
-
-void main(){
-    gl_Position = uProj * vec4(iPos.xyz, 1);
-}

main/src/main.cpp

@@ -1,361 +0,0 @@
-//
-// Created by allem on 2/7/2019.
-//
-
-#include <glad/glad.h>
-#include <GLFW/glfw3.h>
-
-#define GLM_FORCE_SWIZZLE
-#define GLM_ENABLE_EXPERIMENTAL
-
-#include <glm/glm.hpp>
-#include <glm/gtx/string_cast.hpp>
-#include <glm/gtc/matrix_transform.hpp>
-
-#include <cstdlib>
-#include <math.h>
-#include <string>
-#include <vector>
-
-#include "util.h"
-
-#define RES_MAJOR 512
-#define RES_MINOR 8
-
-
-#define Z_RADIUS 10.f
-
-#define PI (float) (M_PI)
-
-struct Unifs {
-    glm::mat4 proj = glm::identity<glm::mat4>();
-    glm::vec4 color = glm::vec4(1);
-};
-
-struct State {
-    GLuint vao_wide{};
-    GLuint vbo_wide{};
-    GLuint ibo_wide{};
-
-    GLuint vao_thin{};
-    GLuint vbo_thin{};
-    GLuint ibo_thin{};
-
-    GLuint ubo{};
-
-    GLuint prog{};
-    const GLuint UNIF_BINDING_POINT = 1;
-
-    std::vector<glm::vec4> verts_wide{};
-    std::vector<unsigned> inds_wide{};
-
-    std::vector<glm::vec4> verts_thin{};
-    std::vector<unsigned> inds_thin{};
-
-    Unifs unifs;
-
-    glm::vec4 hopf(float xi, float nu, float eta) {
-        return glm::vec4(  // todo find parameterization given quaternion
-            cos((nu + xi) / 2) * sin(eta),
-            sin((nu + xi) / 2) * sin(eta),
-            cos((nu - xi) / 2) * cos(eta),
-            sin((nu - xi) / 2) * cos(eta)
-        );
-    }
-
-    glm::vec3 stereo(float xi, float nu, float eta) {
-        auto h = hopf(xi, nu, eta);
-
-        float xw = 0.8;
-        float yw = -0.1;
-        float zw = 2.5;
-        float xy = 0.0;
-        float yz = 0.5;
-        float zx = 0.0;
-
-        auto r = glm::mat4(
-            glm::vec4(cos(xw), 0, 0, sin(xw)),
-            glm::vec4(0, 1, 0, 0),
-            glm::vec4(0, 0, 1, 0),
-            glm::vec4(-sin(xw), 0, 0, cos(xw))
-        ) * glm::mat4(
-            glm::vec4(1, 0, 0, 0),
-            glm::vec4(0, cos(yw), 0, sin(yw)),
-            glm::vec4(0, 0, 1, 0),
-            glm::vec4(0, -sin(yw), 0, cos(yw))
-        ) * glm::mat4(
-            glm::vec4(1, 0, 0, 0),
-            glm::vec4(0, 1, 0, 0),
-            glm::vec4(0, 0, cos(zw), sin(zw)),
-            glm::vec4(0, 0, -sin(zw), cos(zw))
-        ) * glm::mat4(
-            glm::vec4(cos(xy), sin(xy), 0, 0),
-            glm::vec4(-sin(xy), cos(xy), 0, 0),
-            glm::vec4(0, 0, 1, 0),
-            glm::vec4(0, 0, 0, 1)
-        ) * glm::mat4(
-            glm::vec4(1, 0, 0, 0),
-            glm::vec4(0, cos(yz), sin(yz), 0),
-            glm::vec4(0, -sin(yz), cos(yz), 0),
-            glm::vec4(0, 0, 0, 1)
-        ) * glm::mat4(
-            glm::vec4(cos(zx), 0, -sin(zx), 0),
-            glm::vec4(0, 1, 0, 0),
-            glm::vec4(sin(zx), 0, cos(zx), 0),
-            glm::vec4(0, 0, 0, 1)
-        );
-
-        auto rot = r;
-
-        h = rot * h;
-
-        auto s = h.xyz() / (1 - h.w);
-
-        s /= 2;
-
-        return s;
-    }
-
-    void
-    add_ring(std::vector<glm::vec4> &dest_verts, std::vector<unsigned> &dest_inds, float xi, float eta, float rad) {
-        std::vector<glm::vec3> circle;
-        std::vector<glm::vec3> torus;
-        std::vector<unsigned> ind;
-
-        for (unsigned i = 0; i < RES_MAJOR; ++i) {
-            auto nu = 4 * PI * i / RES_MAJOR;
-
-            auto v = stereo(xi, nu, eta);
-
-            circle.push_back(v);
-        }
-
-        auto center = glm::vec3(0);
-        for (auto v : circle) center += v;
-        center /= (float) RES_MAJOR;
-
-        auto A = stereo(xi, 0 * PI, eta);
-        auto B = stereo(xi, 1 * PI, eta);
-        auto normal = glm::normalize(cross(A - center, B - center));
-
-        for (unsigned i = 0; i < circle.size(); ++i) {
-            auto v = circle[i];
-
-            auto b1 = normal;
-            auto b2 = glm::normalize(v - center);
-
-            for (int j = 0; j < RES_MINOR; ++j) {
-                auto theta = 2 * PI * j / RES_MINOR;
-
-                auto p = v + (cos(theta) * b1 + sin(theta) * b2) * rad;
-
-                ind.push_back(i * RES_MINOR + j);
-                ind.push_back((i + 1) * RES_MINOR + j);
-                ind.push_back(i * RES_MINOR + (j + 1) % RES_MINOR);
-
-                ind.push_back((i + 1) * RES_MINOR + j);
-                ind.push_back((i + 1) * RES_MINOR + (j + 1) % RES_MINOR);
-                ind.push_back(i * RES_MINOR + (j + 1) % RES_MINOR);
-
-                torus.push_back(p);
-            }
-        }
-
-        auto offset = (unsigned) dest_verts.size();
-
-        for (auto v : torus)
-            dest_verts.emplace_back(v, 1);
-
-        for (auto i : ind)
-            dest_inds.push_back(offset + i % (RES_MAJOR * RES_MINOR));
-    }
-
-    void updateUnifs() {
-        glBindBuffer(GL_UNIFORM_BUFFER, ubo);
-        util::bufferData(GL_UNIFORM_BUFFER, unifs, GL_STREAM_DRAW);
-        glBindBuffer(GL_UNIFORM_BUFFER, 0);
-    }
-
-    void regen() {
-        printf("generating\n");
-
-        verts_wide.clear();
-        verts_thin.clear();
-        inds_wide.clear();
-        inds_thin.clear();
-
-        const float XI_R = 28;
-        const float ETA_R = 5;
-        const float ETA_BUF = 0;
-
-        for (int i = 0; i < XI_R; ++i) {
-            float xi = 2 * PI * i / XI_R;
-
-            for (int j = 1; j <= ETA_R - 1; ++j) {
-                float eta = ETA_BUF + (PI / 2 - 2 * ETA_BUF) * j / ETA_R;
-
-                add_ring(verts_wide, inds_wide, xi, eta, .015);
-                add_ring(verts_thin, inds_thin, xi, eta, .004);
-            }
-        }
-
-        add_ring(verts_wide, inds_wide, .001f, .0f, .015);
-        add_ring(verts_thin, inds_thin, .001f, .0f, .008);
-
-        add_ring(verts_wide, inds_wide, .001f, PI / 2 - .0f, .015);
-        add_ring(verts_thin, inds_thin, .001f, PI / 2 - .0f, .008);
-
-        glBindBuffer(GL_ARRAY_BUFFER, vbo_thin);
-        util::bufferData(GL_ARRAY_BUFFER, verts_thin, GL_STATIC_DRAW);
-        glBindBuffer(GL_ARRAY_BUFFER, 0);
-
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo_thin);
-        util::bufferData(GL_ELEMENT_ARRAY_BUFFER, inds_thin, GL_STATIC_DRAW);
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
-
-        glBindBuffer(GL_ARRAY_BUFFER, vbo_wide);
-        util::bufferData(GL_ARRAY_BUFFER, verts_wide, GL_STATIC_DRAW);
-        glBindBuffer(GL_ARRAY_BUFFER, 0);
-
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo_wide);
-        util::bufferData(GL_ELEMENT_ARRAY_BUFFER, inds_wide, GL_STATIC_DRAW);
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
-    }
-
-    void init(GLFWwindow *window) {
-        GLuint vs = util::buildShader(GL_VERTEX_SHADER, "vs", {"shaders/main.vert"});
-        GLuint fs = util::buildShader(GL_FRAGMENT_SHADER, "fs", {"shaders/main.frag"});
-        prog = util::buildProgram(false, "prog", {vs, fs});
-
-        glGenBuffers(1, &vbo_wide);
-        glGenBuffers(1, &ibo_wide);
-        glGenBuffers(1, &vbo_thin);
-        glGenBuffers(1, &ibo_thin);
-
-        regen();
-
-        GLint pos = glGetAttribLocation(prog, "iPos");
-
-        glGenBuffers(1, &ubo);
-        glBindBufferBase(GL_UNIFORM_BUFFER, UNIF_BINDING_POINT, ubo);
-        updateUnifs();
-
-        glGenVertexArrays(1, &vao_wide);
-        glBindVertexArray(vao_wide);
-        if (pos >= 0) {
-            glEnableVertexAttribArray(0);
-            glBindBuffer(GL_ARRAY_BUFFER, vbo_wide);
-            glVertexAttribPointer((GLuint) pos, 4, GL_FLOAT, GL_FALSE, sizeof(glm::vec4), nullptr);
-            glBindBuffer(GL_ARRAY_BUFFER, 0);
-        }
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo_wide);
-        glBindVertexArray(0);
-
-        glGenVertexArrays(1, &vao_thin);
-        glBindVertexArray(vao_thin);
-        if (pos >= 0) {
-            glEnableVertexAttribArray(0);
-            glBindBuffer(GL_ARRAY_BUFFER, vbo_thin);
-            glVertexAttribPointer((GLuint) pos, 4, GL_FLOAT, GL_FALSE, sizeof(glm::vec4), nullptr);
-            glBindBuffer(GL_ARRAY_BUFFER, 0);
-        }
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo_thin);
-        glBindVertexArray(0);
-    }
-
-    void render(GLFWwindow *window, float dt, int frame) {
-        int w, h;
-        glfwGetFramebufferSize(window, &w, &h);
-
-        glViewport(0, 0, w, h);
-        glClearColor(1, 1, 1, 1);
-        unifs.color = glm::vec4(0, 0, 0, 1);
-
-        glClear(GL_DEPTH_BUFFER_BIT);
-
-        glBindVertexArray(vao_wide);
-        glEnable(GL_DEPTH_TEST);
-        glEnable(GL_CULL_FACE);
-        glCullFace(GL_BACK);
-        updateUnifs();
-        glUseProgram(prog);
-        glDrawElements(GL_TRIANGLES, (GLsizei) inds_wide.size(), GL_UNSIGNED_INT, 0);
-
-        glClear(GL_COLOR_BUFFER_BIT);
-
-        glBindVertexArray(vao_thin);
-        glEnable(GL_DEPTH_TEST);
-        glEnable(GL_CULL_FACE);
-        glCullFace(GL_BACK);
-        updateUnifs();
-        glUseProgram(prog);
-        glDrawElements(GL_TRIANGLES, (GLsizei) inds_thin.size(), GL_UNSIGNED_INT, 0);
-
-        glfwSwapBuffers(window);
-    }
-
-    void update(GLFWwindow *window, float dt, int frame) {
-        int w, h;
-        glfwGetFramebufferSize(window, &w, &h);
-
-        float ar = (float) w / (float) h;
-        unifs.proj = glm::ortho(-ar, ar, -1.f, 1.f, -Z_RADIUS, Z_RADIUS);
-    }
-
-    void deinit(GLFWwindow *window) {
-
-    }
-
-};
-
-void run(State *state, const std::string &title) {
-    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4);
-    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 4);
-    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
-
-    auto window = glfwCreateWindow(3840, 1249, title.c_str(), nullptr, nullptr);
-    if (!window) {
-        glfwTerminate();
-        exit(EXIT_FAILURE);
-    }
-
-    glfwMakeContextCurrent(window);
-    gladLoadGLLoader((GLADloadproc) glfwGetProcAddress);
-    glfwSwapInterval(0);
-
-    glfwSetWindowUserPointer(window, state);
-    state->init(window);
-
-    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;
-    }
-
-    auto *state = new State();
-    run(state, "Hopf Fibration");
-    delete state;
-
-    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;
+}