5d working

vis/include/mirror.hpp

@@ -1,18 +1,75 @@
 #pragma once
 
 #include <tc/core.hpp>
-#include <glm/glm.hpp>
+
 #include <cmath>
 #include <vector>
 #include <algorithm>
 
-glm::vec4 round(const glm::vec4 &f, int prec) {
-    auto dec = (float) pow(10, prec);
-    auto res = glm::trunc(f * dec + 0.5f) / dec;
-    return res;
+template<unsigned N>
+using vec = std::array<float, N>;
+
+using vec1 = vec<1>;
+using vec2 = vec<2>;
+using vec3 = vec<3>;
+using vec4 = vec<4>;
+using vec5 = vec<5>;
+
+template<class V>
+V operator*(V a, float b) {
+    for (auto &e : a) e *= b;
+    return a;
 }
 
-float dot(int n, const glm::vec4 &a, const glm::vec4 &b) {
+template<class V>
+V operator*(float b, V a) {
+    for (auto &e : a) e *= b;
+    return a;
+}
+
+template<class V>
+V operator/(V a, float b) {
+    for (auto &e : a) e /= b;
+    return a;
+}
+
+template<class V>
+V operator+(V a, V b) {
+    for (int i = 0; i < a.size(); ++i) {
+        a[i] += b[i];
+    }
+    return a;
+}
+
+template<class V>
+V operator-(V a, V b) {
+    for (int i = 0; i < a.size(); ++i) {
+        a[i] -= b[i];
+    }
+    return a;
+}
+
+template<class V>
+void operator-=(V &a, V b) {
+    for (int i = 0; i < a.size(); ++i) {
+        a[i] -= b[i];
+    }
+}
+
+template<class V>
+float length(V a) {
+    float sum = 0;
+    for (auto e : a) sum += e * e;
+    return sqrtf(sum);
+}
+
+template<class V>
+V normalized(V a) {
+    return a / length(a);
+}
+
+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];
@@ -20,15 +77,17 @@ float dot(int n, const glm::vec4 &a, const glm::vec4 &b) {
     return sum;
 }
 
-float dot(int n, const std::vector<float> &a, const std::vector<float> &b) {
+template<class V>
+float dot(const V &a, const V &b) {
     float sum = 0;
-    for (int i = 0; i < n; ++i) {
+    for (int i = 0; i < a.size(); ++i) {
         sum += a[i] * b[i];
     }
     return sum;
 }
 
-std::vector<glm::vec4> mirror(const tc::Group &group) {
+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) {
@@ -51,33 +110,41 @@ std::vector<glm::vec4> mirror(const tc::Group &group) {
         mirrors.push_back(vp);
     }
 
-    std::vector<glm::vec4> res;
-    for (const auto &vec : mirrors) {
-        glm::vec4 rvec{};
+    std::vector<vec<N>> res;
+    for (const auto &v : mirrors) {
+        vec<N> rv{};
 
         // ortho proj
-        for (int i = 0; i < std::min(vec.size(), (size_t) 4); ++i) {
-            rvec[i] = vec[i];
+        for (int i = 0; i < std::min(v.size(), (size_t) N); ++i) {
+            rv[i] = v[i];
         }
 
-        res.push_back(rvec);
+        res.push_back(rv);
     }
     return res;
 }
 
-glm::vec4 project(const glm::vec4 &vec, const glm::vec4 &target) {
-    return glm::dot(vec, target) / glm::dot(target, target) * target;
+template<unsigned N>
+vec<N> stereo(vec<N + 1> v) {
+    vec<N> r;
+    for (int i = 0; i < N; ++i) {
+        r[i] = v[i] / (1-v[N]);
+    }
+    return r;
 }
 
-glm::vec4 reflect(const glm::vec4 &vec, const glm::vec4 axis) {
-    return vec - 2.f * project(vec, axis);
+template<class V>
+V project(const V &vec, const V &target) {
+    return dot(vec, target) / dot(target, target) * target;
 }
 
-glm::vec4 reflect_scaled(const  glm::vec4 &vec, const glm::vec4 axis) {
-    return vec - 2.f * glm::length(axis) * project(vec, axis);
+template<class V>
+V reflect(const V &a, const V &axis) {
+    return a - 2.f * project(a, axis);
 }
 
-glm::vec4 gram_schmidt_last(std::vector<glm::vec4> vecs) {
+template<class V>
+V gram_schmidt_last(std::vector<V> vecs) {
     int N = vecs.size();
     for (int i = 0; i < N; ++i) {
         for (int j = 0; j < i; ++j) {
@@ -85,22 +152,24 @@ glm::vec4 gram_schmidt_last(std::vector<glm::vec4> vecs) {
         }
     }
 
-    return glm::normalize(vecs[N - 1]);
+    return normalized(vecs[N - 1]);
 }
 
-glm::vec4 barycentric(std::vector<glm::vec4> basis, std::vector<float> coords) {
-    glm::vec4 res{};
+template<class V>
+V barycentric(std::vector<V> basis, std::vector<float> coords) {
+    V res{};
 
     int N = std::min(basis.size(), coords.size());
     for (int i = 0; i < N; ++i) {
-        res += basis[i] * coords[i];
+        res = res + (basis[i] * coords[i]);
     }
-    return glm::normalize(res);
+    return normalized(res);
 }
 
-std::vector<glm::vec4> plane_intersections(std::vector<glm::vec4> normals) {
+template<class V>
+std::vector<V> plane_intersections(std::vector<V> normals) {
     int N = normals.size();
-    std::vector<glm::vec4> results(N);
+    std::vector<V> results(N);
 
     for (int i = 0; i < N; ++i) {
         std::rotate(normals.begin(), normals.begin() + 1, normals.end());