ENH: Replace std::array<float, N> with Eigen

vis/include/mirror.hpp

@@ -6,76 +6,7 @@
 #include <vector>
 #include <algorithm>
 
-#include <Eigen/Eigen>
-
-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, const float &b) {
-    for (auto &e : a) e *= b;
-    return a;
-}
-
-template<class V>
-V operator*(const float &b, V a) {
-    for (auto &e : a) e *= b;
-    return a;
-}
-
-template<class V>
-V operator/(V a, const float &b) {
-    for (auto &e : a) e /= b;
-    return a;
-}
-
-template<class V>
-V operator+(const 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, const V &b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] -= b[i];
-    }
-    return a;
-}
-
-template<class V>
-void operator-=(V &a, const V &b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] -= b[i];
-    }
-}
-
-template<class V>
-void operator+=(V &a, const V &b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] += b[i];
-    }
-}
-
-template<class V>
-float length(const V &a) {
-    float sum = 0;
-    for (const auto &e : a) sum += e * e;
-    return sqrtf(sum);
-}
-
-template<class V>
-V normalized(const V &a) {
-    return a / length(a);
-}
+#include <geometry.hpp>
 
 template<class V>
 float dot(int n, const V &a, const V &b) {
@@ -86,15 +17,6 @@ float dot(int n, const V &a, const V &b) {
     return sum;
 }
 
-template<class V>
-float dot(const V &a, const V &b) {
-    float sum = 0;
-    for (int i = 0; i < a.size(); ++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;
@@ -121,7 +43,7 @@ std::vector<vec<N>> mirror(const tc::Group<> &group) {
 
     std::vector<vec<N>> res;
     for (const auto &v : mirrors) {
-        vec<N> rv{};
+        vec<N> rv = vec<N>::Zero();
 
         // ortho proj
         for (int i = 0; i < std::min(v.size(), (size_t) N); ++i) {
@@ -142,9 +64,18 @@ vec<N> stereo(const vec<N + 1> &v) {
     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 dot(vec, target) / dot(target, target) * target;
+    return vec.dot(target) / target.dot(target) * target;
 }
 
 template<class V>
@@ -160,14 +91,14 @@ V gram_schmidt_last(std::vector<V> vecs) {
         }
     }
 
-    return normalized(vecs[vecs.size() - 1]);
+    return vecs[vecs.size() - 1].normalized();
 }
 
 template<class V, class C>
 V barycentric(const std::vector<V> &basis, const C &coords) {
-    V res{};
+    V res = V::Zero();
 
-    int N = std::min(basis.size(), coords.size());
+    int N = std::min((int) basis.size(), (int) coords.rows());
     for (int i = 0; i < N; ++i) {
         res += basis[i] * coords[i];
     }
@@ -186,29 +117,12 @@ std::vector<V> plane_intersections(std::vector<V> normals) {
     return results;
 }
 
-Eigen::Matrix4f utilRotate(const int u, const int v, const float theta) {
-    Eigen::Matrix4f res;
-    res.setIdentity();
+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;
 }
-
-Eigen::Matrix4f ortho(
-    float l,
-    float r,
-    float b,
-    float t,
-    float n,
-    float f
-) {
-    Eigen::Matrix4f res;
-    res <<
-        2 / (r - l), 0, 0, -(r + l) / (r - l),
-        0, 2 / (t - b), 0, -(t + b) / (t - b),
-        0, 0, -2 / (f - n), -(f + n) / (f - n),
-        0, 0, 0, 1;
-    return res;
-}