diff --git a/bsputil/bsputil.cc b/bsputil/bsputil.cc
index 5ea3a8c3..9b378ec2 100644
--- a/bsputil/bsputil.cc
+++ b/bsputil/bsputil.cc
@@ -444,7 +444,7 @@ static void CompareBSPFiles(const mbsp_t &refBsp, const mbsp_t &bsp)
     }
 }
 
-static void FindFaces(const mbsp_t *bsp, const vec3_t &pos, const vec3_t &normal)
+static void FindFaces(const mbsp_t *bsp, const qvec3d &pos, const qvec3d &normal)
 {
     for (int i = 0; i < bsp->dmodels.size(); ++i) {
         const dmodelh2_t *model = &bsp->dmodels[i];
@@ -567,8 +567,8 @@ int main(int argc, char **argv)
             }
 
             try {
-                const vec3_t pos = {std::stof(argv[i + 1]), std::stof(argv[i + 2]), std::stof(argv[i + 3])};
-                const vec3_t normal = {std::stof(argv[i + 4]), std::stof(argv[i + 5]), std::stof(argv[i + 6])};
+                const qvec3d pos {std::stof(argv[i + 1]), std::stof(argv[i + 2]), std::stof(argv[i + 3])};
+                const qvec3d normal {std::stof(argv[i + 4]), std::stof(argv[i + 5]), std::stof(argv[i + 6])};
                 FindFaces(&bsp, pos, normal);
             }
             catch (const std::exception &) {
diff --git a/bsputil/decompile.cpp b/bsputil/decompile.cpp
index e24028a5..ccd46943 100644
--- a/bsputil/decompile.cpp
+++ b/bsputil/decompile.cpp
@@ -145,9 +145,7 @@ std::vector<decomp_plane_t> RemoveRedundantPlanes(const std::vector<decomp_plane
 
     for (const decomp_plane_t &plane : planes) {
         // outward-facing plane
-        vec3_t normal;
-        VectorCopy(plane.normal, normal);
-        std::optional<winding_t> winding = winding_t::from_plane(normal, plane.distance, 10e6);
+        std::optional<winding_t> winding = winding_t::from_plane(plane.normal, plane.distance, 10e6);
 
         // clip `winding` by all of the other planes, flipped
         for (const decomp_plane_t &plane2 : planes) {
@@ -302,10 +300,7 @@ public:
      */
     std::pair<decomp_brush_face_t, decomp_brush_face_t> clipToPlane(const qvec3d &normal, double distance) const
     {
-        vec3_t pnormal;
-        VectorCopy(normal, pnormal);
-
-        auto clipped = winding->clip(pnormal, (float)distance);
+        auto clipped = winding->clip(normal, (float)distance);
 
         // front or back may be null (if fully clipped).
         // these constructors take ownership of the winding.
@@ -315,8 +310,7 @@ public:
 
     qvec3d normal() const
     {
-        plane_t plane = winding->plane();
-        return plane.normal;
+        return winding->plane().normal;
     }
 };
 
diff --git a/common/bsputils.cc b/common/bsputils.cc
index cbc0f9a8..9932ba70 100644
--- a/common/bsputils.cc
+++ b/common/bsputils.cc
@@ -126,17 +126,15 @@ qplane3d Face_Plane(const mbsp_t *bsp, const mface_t *f)
     Q_assert(f->planenum >= 0 && f->planenum < bsp->dplanes.size());
     const dplane_t *dplane = &bsp->dplanes[f->planenum];
 
-    vec3_t planeNormal;
-    VectorCopy(dplane->normal, planeNormal); // convert from float->double if needed
-
     plane_t result;
     if (f->side) {
-        VectorSubtract(vec3_origin, planeNormal, result.normal);
+        result.normal = -dplane->normal;
         result.dist = -dplane->dist;
     } else {
-        VectorCopy(planeNormal, result.normal);
+        result.normal = dplane->normal;
         result.dist = dplane->dist;
     }
+
     return {result.normal, result.dist};
 }
 
@@ -276,9 +274,7 @@ vec_t Plane_Dist(const qvec3d &point, const dplane_t *plane)
         case PLANE_Y: return point[1] - plane->dist;
         case PLANE_Z: return point[2] - plane->dist;
         default: {
-            vec3_t planeNormal;
-            VectorCopy(plane->normal, planeNormal); // convert from float->double if needed
-            return DotProduct(point, planeNormal) - plane->dist;
+            return DotProduct(point, qvec3d(plane->normal)) - plane->dist;
         }
     }
 }
@@ -388,12 +384,10 @@ plane_t *Face_AllocInwardFacingEdgePlanes(const mbsp_t *bsp, const mface_t *face
     for (int i = 0; i < face->numedges; i++) {
         plane_t *dest = &out[i];
 
-        const qvec3d v0 = GetSurfaceVertexPoint(bsp, face, i);
-        const qvec3d v1 = GetSurfaceVertexPoint(bsp, face, (i + 1) % face->numedges);
+        const qvec3d &v0 = GetSurfaceVertexPoint(bsp, face, i);
+        const qvec3d &v1 = GetSurfaceVertexPoint(bsp, face, (i + 1) % face->numedges);
 
-        vec3_t edgevec;
-        VectorSubtract(v1, v0, edgevec);
-        VectorNormalize(edgevec);
+        qvec3d edgevec = qv::normalize(v1 - v0);
 
         CrossProduct(edgevec, faceplane.normal(), dest->normal);
         dest->dist = DotProduct(dest->normal, v0);
diff --git a/include/common/qvec.hh b/include/common/qvec.hh
index ed72667f..8c4d140f 100644
--- a/include/common/qvec.hh
+++ b/include/common/qvec.hh
@@ -130,12 +130,12 @@ public:
     [[nodiscard]] constexpr size_t size() const { return N; }
 
     // Sort support
-    [[nodiscard]] constexpr bool operator<(const qvec<T, N> &other) const { return v < other.v; }
-    [[nodiscard]] constexpr bool operator<=(const qvec<T, N> &other) const { return v <= other.v; }
-    [[nodiscard]] constexpr bool operator>(const qvec<T, N> &other) const { return v > other.v; }
-    [[nodiscard]] constexpr bool operator>=(const qvec<T, N> &other) const { return v >= other.v; }
-    [[nodiscard]] constexpr bool operator==(const qvec<T, N> &other) const { return v == other.v; }
-    [[nodiscard]] constexpr bool operator!=(const qvec<T, N> &other) const { return v != other.v; }
+    [[nodiscard]] constexpr bool operator<(const qvec &other) const { return v < other.v; }
+    [[nodiscard]] constexpr bool operator<=(const qvec &other) const { return v <= other.v; }
+    [[nodiscard]] constexpr bool operator>(const qvec &other) const { return v > other.v; }
+    [[nodiscard]] constexpr bool operator>=(const qvec &other) const { return v >= other.v; }
+    [[nodiscard]] constexpr bool operator==(const qvec &other) const { return v == other.v; }
+    [[nodiscard]] constexpr bool operator!=(const qvec &other) const { return v != other.v; }
 
     [[nodiscard]] constexpr const T &at(const size_t idx) const
     {
@@ -153,63 +153,79 @@ public:
     [[nodiscard]] constexpr T &operator[](const size_t idx) { return at(idx); }
 
     template<typename F>
-    constexpr void operator+=(const qvec<F, N> &other)
+    constexpr qvec operator+=(const qvec<F, N> &other)
     {
         for (size_t i = 0; i < N; i++)
             v[i] += other.v[i];
+        return *this;
     }
     template<typename F>
-    constexpr void operator-=(const qvec<F, N> &other)
+    constexpr qvec operator-=(const qvec<F, N> &other)
     {
         for (size_t i = 0; i < N; i++)
             v[i] -= other.v[i];
+        return *this;
     }
-    constexpr void operator*=(const T &scale)
+    constexpr qvec operator*=(const T &scale)
     {
         for (size_t i = 0; i < N; i++)
             v[i] *= scale;
+        return *this;
     }
-    constexpr void operator/=(const T &scale)
+    constexpr qvec operator/=(const T &scale)
     {
         for (size_t i = 0; i < N; i++)
             v[i] /= scale;
+        return *this;
+    }
+    constexpr qvec &operator*=(const qvec &scale)
+    {
+        for (size_t i = 0; i < N; i++)
+            v[i] *= scale[i];
+        return *this;
+    }
+    constexpr qvec operator/=(const qvec &scale)
+    {
+        for (size_t i = 0; i < N; i++)
+            v[i] /= scale[i];
+        return *this;
     }
 
     template<typename F>
-    [[nodiscard]] constexpr qvec<T, N> operator+(const qvec<F, N> &other) const
+    [[nodiscard]] constexpr qvec operator+(const qvec<F, N> &other) const
     {
-        qvec<T, N> res(*this);
-        res += other;
-        return res;
+        return qvec(*this) += other;
     }
 
     template<typename F>
     [[nodiscard]] constexpr qvec<T, N> operator-(const qvec<F, N> &other) const
     {
-        qvec<T, N> res(*this);
-        res -= other;
-        return res;
+        return qvec(*this) -= other;
     }
 
     [[nodiscard]] constexpr qvec<T, N> operator*(const T &scale) const
     {
-        qvec<T, N> res(*this);
-        res *= scale;
-        return res;
+        return qvec(*this) *= scale;
     }
 
     [[nodiscard]] constexpr qvec<T, N> operator/(const T &scale) const
     {
-        qvec<T, N> res(*this);
-        res /= scale;
-        return res;
+        return qvec(*this) /= scale;
+    }
+    
+    [[nodiscard]] constexpr qvec<T, N> operator*(const qvec<T, N> &scale) const
+    {
+        return qvec(*this) *= scale;
+    }
+
+    [[nodiscard]] constexpr qvec<T, N> operator/(const qvec<T, N> &scale) const
+    {
+        return qvec(*this) /= scale;
     }
 
     [[nodiscard]] constexpr qvec<T, N> operator-() const
     {
-        qvec<T, N> res(*this);
-        res *= -1;
-        return res;
+        return qvec(*this) *= -1;
     }
 
     [[nodiscard]] constexpr qvec<T, 3> xyz() const
@@ -230,7 +246,7 @@ public:
     constexpr auto cend() const { return v.cend(); }
 
     // for Google Test
-    friend std::ostream& operator<<(std::ostream& os, const qvec<T, N>& v) {
+    friend std::ostream &operator<<(std::ostream &os, const qvec &v) {
         return os << fmt::format("{}", v);
     }
 };
@@ -501,6 +517,26 @@ template<typename T>
 {
     return (t0 * bary[0]) + (t1 * bary[1]) + (t2 * bary[2]);
 }
+
+// Snap vector to nearest axial component
+template<typename T>
+[[nodiscard]] qvec<T, 3> Snap(const qvec<T, 3> &normal, const T &epsilon = NORMAL_EPSILON)
+{
+    qvec<T, 3> vec { };
+
+    for (size_t i = 0; i < 3; i++) {
+        if (fabs(normal[i] - 1) < epsilon) {
+            vec[i] = 1;
+            break;
+        }
+        if (fabs(normal[i] - -1) < epsilon) {
+            vec[i] = -1;
+            break;
+        }
+    }
+    
+    return vec;
+}
 }; // namespace qv
 
 using qvec2f = qvec<float, 2>;
@@ -532,6 +568,16 @@ public:
     {
     }
 
+    auto as_tuple() const { return std::tie(m_normal, m_dist); }
+
+    // Sort support
+    [[nodiscard]] constexpr bool operator<(const qplane3 &other) const { return as_tuple() < other.as_tuple(); }
+    [[nodiscard]] constexpr bool operator<=(const qplane3 &other) const { return as_tuple() <= other.as_tuple(); }
+    [[nodiscard]] constexpr bool operator>(const qplane3 &other) const { return as_tuple() > other.as_tuple(); }
+    [[nodiscard]] constexpr bool operator>=(const qplane3 &other) const { return as_tuple() >= other.as_tuple(); }
+    [[nodiscard]] constexpr bool operator==(const qplane3 &other) const { return as_tuple() == other.as_tuple(); }
+    [[nodiscard]] constexpr bool operator!=(const qplane3 &other) const { return as_tuple() != other.as_tuple(); }
+
     [[nodiscard]] inline T distAbove(const qvec<T, 3> &pt) const { return qv::dot(pt, m_normal) - m_dist; }
     [[nodiscard]] constexpr const qvec<T, 3> &normal() const { return m_normal; }
     [[nodiscard]] constexpr const T dist() const { return m_dist; }
@@ -546,16 +592,16 @@ using qplane3f = qplane3<float>;
 using qplane3d = qplane3<double>;
 
 /**
- * M row, N column matrix.
+ * Row x Col matrix of T.
  */
-template<size_t M, size_t N, class T>
+template<class T, size_t NRow, size_t NCol>
 class qmat
 {
 public:
     /**
      * Column-major order. [ (row0,col0), (row1,col0), .. ]
      */
-    std::array<T, M * N> m_values;
+    std::array<T, NRow * NCol> m_values;
 
 public:
     /**
@@ -563,9 +609,9 @@ public:
      */
     constexpr qmat() : m_values({})
     {
-        if constexpr (M == N) {
+        if constexpr (NRow == NCol) {
             // identity matrix
-            for (size_t i = 0; i < N; i++) {
+            for (size_t i = 0; i < NCol; i++) {
                 this->at(i, i) = 1;
             }
         }
@@ -577,63 +623,73 @@ public:
     inline qmat(const T &val) { m_values.fill(val); }
 
     // copy constructor
-    constexpr qmat(const qmat<M, N, T> &other) : m_values(other.m_values) { }
+    constexpr qmat(const qmat &other) : m_values(other.m_values) { }
 
     /**
      * Casting from another matrix type of the same size
      */
     template<class T2>
-    constexpr qmat(const qmat<M, N, T2> &other)
+    constexpr qmat(const qmat<T2, NRow, NCol> &other)
     {
-        for (size_t i = 0; i < M * N; i++)
+        for (size_t i = 0; i < NRow * NCol; i++)
             this->m_values[i] = static_cast<T>(other.m_values[i]);
     }
 
     // initializer list, column-major order
     constexpr qmat(std::initializer_list<T> list)
     {
-        assert(list.size() == M * N);
+        assert(list.size() == NRow * NCol);
         std::copy(list.begin(), list.end(), m_values.begin());
     }
-
-    constexpr bool operator==(const qmat<M, N, T> &other) const { return m_values == other.m_values; }
+    
+    // Sort support
+    [[nodiscard]] constexpr bool operator<(const qmat &other) const { return m_values < other.m_values; }
+    [[nodiscard]] constexpr bool operator<=(const qmat &other) const { return m_values <= other.m_values; }
+    [[nodiscard]] constexpr bool operator>(const qmat &other) const { return m_values > other.m_values; }
+    [[nodiscard]] constexpr bool operator>=(const qmat &other) const { return m_values >= other.m_values; }
+    [[nodiscard]] constexpr bool operator==(const qmat &other) const { return m_values == other.m_values; }
+    [[nodiscard]] constexpr bool operator!=(const qmat &other) const { return m_values != other.m_values; }
 
     // access to elements
 
     [[nodiscard]] constexpr T &at(size_t row, size_t col)
     {
-        assert(row >= 0 && row < M);
-        assert(col >= 0 && col < N);
-        return m_values[col * M + row];
+        assert(row >= 0 && row < NRow);
+        assert(col >= 0 && col < NCol);
+        return m_values[col * NRow + row];
     }
 
     [[nodiscard]] constexpr T at(size_t row, size_t col) const
     {
-        assert(row >= 0 && row < M);
-        assert(col >= 0 && col < N);
-        return m_values[col * M + row];
+        assert(row >= 0 && row < NRow);
+        assert(col >= 0 && col < NCol);
+        return m_values[col * NRow + row];
     }
 
-    // hacky accessor for mat[col][row] access
-    [[nodiscard]] constexpr const T *operator[](size_t col) const
+    // access row
+    [[nodiscard]] constexpr qvec<T, NCol> row(size_t row) const
     {
-        assert(col >= 0 && col < N);
-        return &m_values[col * M];
+        assert(row >= 0 && row < NRow);
+        qvec<T, NCol> v;
+        for (size_t i = 0; i < NCol; i++) {
+            v[i] = at(row, i);
+        }
+        return v;
     }
 
-    [[nodiscard]] constexpr T *operator[](size_t col)
+    [[nodiscard]] constexpr const qvec<T, NRow> &col(size_t col) const
     {
-        assert(col >= 0 && col < N);
-        return &m_values[col * M];
+        assert(col >= 0 && col < NCol);
+        return reinterpret_cast<const qvec<T, NRow> &>(m_values[col * NRow]);
     }
 
     // multiplication by a vector
 
-    [[nodiscard]] constexpr qvec<T, M> operator*(const qvec<T, N> &vec) const
+    [[nodiscard]] constexpr qvec<T, NRow> operator*(const qvec<T, NCol> &vec) const
     {
-        qvec<T, M> res{};
-        for (size_t i = 0; i < M; i++) { // for each row
-            for (size_t j = 0; j < N; j++) { // for each col
+        qvec<T, NRow> res{};
+        for (size_t i = 0; i < NRow; i++) { // for each row
+            for (size_t j = 0; j < NCol; j++) { // for each col
                 res[i] += this->at(i, j) * vec[j];
             }
         }
@@ -642,14 +698,14 @@ public:
 
     // multiplication by a matrix
 
-    template<size_t P>
-    [[nodiscard]] constexpr qmat<M, P, T> operator*(const qmat<N, P, T> &other) const
+    template<size_t PCol>
+    [[nodiscard]] constexpr qmat<T, NRow, PCol> operator*(const qmat<T, NCol, PCol> &other) const
     {
-        qmat<M, P, T> res;
-        for (size_t i = 0; i < M; i++) {
-            for (size_t j = 0; j < P; j++) {
+        qmat<T, NRow, PCol> res;
+        for (size_t i = 0; i < NRow; i++) {
+            for (size_t j = 0; j < PCol; j++) {
                 T val = 0;
-                for (size_t k = 0; k < N; k++) {
+                for (size_t k = 0; k < NCol; k++) {
                     val += this->at(i, k) * other.at(k, j);
                 }
                 res.at(i, j) = val;
@@ -660,23 +716,26 @@ public:
 
     // multiplication by a scalar
 
-    [[nodiscard]] constexpr qmat<M, N, T> operator*(const T scalar) const
+    [[nodiscard]] constexpr qmat operator*(const T &scalar) const
     {
-        qmat<M, N, T> res(*this);
-        for (size_t i = 0; i < M * N; i++) {
+        qmat res(*this);
+        for (size_t i = 0; i < NRow * NCol; i++) {
             res.m_values[i] *= scalar;
         }
         return res;
     }
+
+    // stream support
+    auto stream_data() { return std::tie(m_values); }
 };
 
-using qmat2x2f = qmat<2, 2, float>;
-using qmat3x3f = qmat<3, 3, float>;
-using qmat4x4f = qmat<4, 4, float>;
+using qmat2x2f = qmat<float, 2, 2>;
+using qmat3x3f = qmat<float, 3, 3>;
+using qmat4x4f = qmat<float, 4, 4>;
 
-using qmat2x2d = qmat<2, 2, double>;
-using qmat3x3d = qmat<3, 3, double>;
-using qmat4x4d = qmat<4, 4, double>;
+using qmat2x2d = qmat<double, 2, 2>;
+using qmat3x3d = qmat<double, 3, 3>;
+using qmat4x4d = qmat<double, 4, 4>;
 
 namespace qv
 {
diff --git a/include/qbsp/brush.hh b/include/qbsp/brush.hh
index 9ecd8146..2b103d0b 100644
--- a/include/qbsp/brush.hh
+++ b/include/qbsp/brush.hh
@@ -52,13 +52,13 @@ brush_t *LoadBrush(const mapentity_t *src, const mapbrush_t *mapbrush, const con
     const qvec3d &rotate_offset, const rotation_t rottype, const int hullnum);
 void FreeBrushes(mapentity_t *ent);
 
-int FindPlane(const vec3_t normal, const vec_t dist, int *side);
+int FindPlane(const qvec3d &normal, const vec_t dist, int *side);
 bool PlaneEqual(const plane_t *p1, const plane_t *p2);
 bool PlaneInvEqual(const plane_t *p1, const plane_t *p2);
 
 bool BoundBrush(brush_t *brush);
 vec_t BrushVolume(const brush_t *brush);
-int BrushMostlyOnSide(const brush_t *brush, const vec3_t normal, vec_t dist);
+int BrushMostlyOnSide(const brush_t *brush, const qvec3d &normal, vec_t dist);
 void SplitBrush(const brush_t *brush, int planenum, int planeside, brush_t **front, brush_t **back);
 
 void FilterStructuralBrushesIntoTree(const mapentity_t *e, node_t *headnode);
diff --git a/include/qbsp/qbsp.hh b/include/qbsp/qbsp.hh
index 34c43713..ac80e4c2 100644
--- a/include/qbsp/qbsp.hh
+++ b/include/qbsp/qbsp.hh
@@ -111,6 +111,7 @@ enum
 #include <common/mathlib.hh>
 #include <qbsp/winding.hh>
 
+
 struct mtexinfo_t
 {
     texvecf vecs; /* [s/t][xyz offset] */
@@ -124,7 +125,6 @@ struct mtexinfo_t
     constexpr bool operator<(const mtexinfo_t &other) const { return as_tuple() < other.as_tuple(); }
 
     constexpr bool operator>(const mtexinfo_t &other) const { return as_tuple() > other.as_tuple(); }
-
     constexpr qvec2d uvs(const qvec3d &pos, const int32_t &width, const int32_t &height) const
     {
         return {(pos[0] * vecs[0][0] + pos[1] * vecs[0][1] + pos[2] * vecs[0][2] + vecs[0][3]) / width,
diff --git a/light/light.cc b/light/light.cc
index a97afc34..62f6f5bb 100644
--- a/light/light.cc
+++ b/light/light.cc
@@ -904,9 +904,9 @@ static inline void WriteNormals(const mbsp_t &bsp, bspdata_t &bspdata)
     for (auto &face : bsp.dfaces) {
         auto &cache = FaceCacheForFNum(&face - bsp.dfaces.data());
         for (auto &normals : cache.normals()) {
-            unique_normals.insert(normals.normal);
-            unique_normals.insert(normals.tangent);
-            unique_normals.insert(normals.bitangent);
+            unique_normals.insert(qv::Snap(normals.normal));
+            unique_normals.insert(qv::Snap(normals.tangent));
+            unique_normals.insert(qv::Snap(normals.bitangent));
             num_normals += 3;
         }
     }
@@ -929,9 +929,9 @@ static inline void WriteNormals(const mbsp_t &bsp, bspdata_t &bspdata)
         auto &cache = FaceCacheForFNum(&face - bsp.dfaces.data());
 
         for (auto &n : cache.normals()) {
-            stream <= numeric_cast<uint32_t>(mapped_normals[n.normal]);
-            stream <= numeric_cast<uint32_t>(mapped_normals[n.tangent]);
-            stream <= numeric_cast<uint32_t>(mapped_normals[n.bitangent]);
+            stream <= numeric_cast<uint32_t>(mapped_normals[qv::Snap(n.normal)]);
+            stream <= numeric_cast<uint32_t>(mapped_normals[qv::Snap(n.tangent)]);
+            stream <= numeric_cast<uint32_t>(mapped_normals[qv::Snap(n.bitangent)]);
         }
     }
 
diff --git a/qbsp/brush.cc b/qbsp/brush.cc
index 19ba4fbd..71892f18 100644
--- a/qbsp/brush.cc
+++ b/qbsp/brush.cc
@@ -218,26 +218,23 @@ static void PlaneHash_Add(const qbsp_plane_t *p, int index)
  * NewPlane
  * - Returns a global plane number and the side that will be the front
  */
-static int NewPlane(const vec3_t normal, const vec_t dist, int *side)
+static int NewPlane(const qbsp_plane_t &plane, int *side)
 {
-    vec_t len;
+    vec_t len = VectorLength(plane.normal);
 
-    len = VectorLength(normal);
     if (len < 1 - ON_EPSILON || len > 1 + ON_EPSILON)
         FError("invalid normal (vector length {:.4})", len);
 
     size_t index = map.planes.size();
-    qbsp_plane_t &plane = map.planes.emplace_back();
-    VectorCopy(normal, plane.normal);
-    plane.dist = dist;
+    qbsp_plane_t &added_plane = map.planes.emplace_back(plane);
 
-    int32_t out_side = NormalizePlane(&plane, side != nullptr);
+    int32_t out_side = NormalizePlane(&added_plane, side != nullptr);
 
     if (side) {
         *side = out_side;
     }
 
-    PlaneHash_Add(&plane, index);
+    PlaneHash_Add(&added_plane, index);
     return index;
 }
 
@@ -246,10 +243,10 @@ static int NewPlane(const vec3_t normal, const vec_t dist, int *side)
  * - Returns a global plane number and the side that will be the front
  * - if `side` is null, only an exact match will be fetched.
  */
-int FindPlane(const vec3_t normal, const vec_t dist, int *side)
+int FindPlane(const qvec3d &normal, const vec_t dist, int *side)
 {
-    qbsp_plane_t plane = {0};
-    VectorCopy(normal, plane.normal);
+    qbsp_plane_t plane {};
+    plane.normal = normal;
     plane.dist = dist;
 
     for (int i : map.planehash[plane_hash_fn(&plane)]) {
@@ -264,7 +261,7 @@ int FindPlane(const vec3_t normal, const vec_t dist, int *side)
             return i;
         }
     }
-    return NewPlane(&plane.normal[0], plane.dist, side);
+    return NewPlane(plane, side);
 }
 
 /*
@@ -423,7 +420,7 @@ static face_t *CreateBrushFaces(const mapentity_t *src, hullbrush_t *hullbrush,
         plane.dist = DotProduct(plane.normal, point);
 
         f->texinfo = hullnum > 0 ? 0 : mapface->texinfo;
-        f->planenum = FindPlane(&plane.normal[0], plane.dist, &f->planeside);
+        f->planenum = FindPlane(plane.normal, plane.dist, &f->planeside);
         f->next = facelist;
         facelist = f;
         CheckFace(f, mapface);
@@ -1328,7 +1325,7 @@ BrushMostlyOnSide
 from q3map
 ==================
 */
-int BrushMostlyOnSide(const brush_t *brush, const vec3_t planenormal, vec_t planedist)
+int BrushMostlyOnSide(const brush_t *brush, const qvec3d &planenormal, vec_t planedist)
 {
     vec_t max;
     int side;
@@ -1523,7 +1520,7 @@ void SplitBrush(const brush_t *brush, int planenum, int planeside, brush_t **fro
     if (!w || WindingIsTiny(*w)) { // the brush isn't really split
         int side;
 
-        side = BrushMostlyOnSide(brush, &plane.normal[0], plane.dist);
+        side = BrushMostlyOnSide(brush, plane.normal, plane.dist);
         if (side == SIDE_FRONT)
             *front = CopyBrush(brush);
         if (side == SIDE_BACK)
diff --git a/qbsp/map.cc b/qbsp/map.cc
index 1ffd6b04..ff688155 100644
--- a/qbsp/map.cc
+++ b/qbsp/map.cc
@@ -835,8 +835,8 @@ static texdef_quake_ed_noshift_t Reverse_QuakeEd(qmat2x2f M, const qbsp_plane_t
 {
     // Check for shear, because we might tweak M to remove it
     {
-        qvec2f Xvec = qvec2f(M[0][0], M[1][0]);
-        qvec2f Yvec = qvec2f(M[0][1], M[1][1]);
+        qvec2f Xvec = M.row(0);
+        qvec2f Yvec = M.row(1);
         double cosAngle = qv::dot(qv::normalize(Xvec), qv::normalize(Yvec));
 
         // const double oldXscale = sqrt(pow(M[0][0], 2.0) + pow(M[1][0], 2.0));
@@ -878,17 +878,17 @@ static texdef_quake_ed_noshift_t Reverse_QuakeEd(qmat2x2f M, const qbsp_plane_t
             }
 
             // update M
-            M[0][0] = Xvec[0];
-            M[1][0] = Xvec[1];
+            M.at(0, 0) = Xvec[0];
+            M.at(0, 1) = Xvec[1];
 
-            M[0][1] = Yvec[0];
-            M[1][1] = Yvec[1];
+            M.at(1, 0) = Yvec[0];
+            M.at(1, 1) = Yvec[1];
         }
     }
 
     // extract abs(scale)
-    const double absXscale = sqrt(pow(M[0][0], 2.0) + pow(M[1][0], 2.0));
-    const double absYscale = sqrt(pow(M[0][1], 2.0) + pow(M[1][1], 2.0));
+    const double absXscale = sqrt(pow(M.at(0, 0), 2.0) + pow(M.at(0, 1), 2.0));
+    const double absYscale = sqrt(pow(M.at(1, 0), 2.0) + pow(M.at(1, 1), 2.0));
     const qmat2x2f applyAbsScaleM{static_cast<float>(absXscale), // col0
         0,
         0, // col1
@@ -933,8 +933,8 @@ static texdef_quake_ed_noshift_t Reverse_QuakeEd(qmat2x2f M, const qbsp_plane_t
             const qmat2x2f applyAngleGuessM = rotation2x2_deg(angleGuess);
             const qmat2x2f Mguess = applyGuessedFlipM * applyAbsScaleM * applyAngleGuessM * axisFlipsM;
 
-            if (fabs(M[0][0] - Mguess[0][0]) < 0.001 && fabs(M[0][1] - Mguess[0][1]) < 0.001 &&
-                fabs(M[1][0] - Mguess[1][0]) < 0.001 && fabs(M[1][1] - Mguess[1][1]) < 0.001) {
+            if (fabs(M.at(0, 0) - Mguess.at(0, 0)) < 0.001 && fabs(M.at(1, 0) - Mguess.at(1, 0)) < 0.001 &&
+                fabs(M.at(0, 1) - Mguess.at(0, 1)) < 0.001 && fabs(M.at(1, 1) - Mguess.at(1, 1)) < 0.001) {
 
                 texdef_quake_ed_noshift_t reversed;
                 reversed.rotate = angleGuess;
@@ -1001,7 +1001,6 @@ static void SetTexinfo_QuakeEd_New(
     }
 
     // copy M into the output vectors
-
     for (int i = 0; i < 2; i++) {
         for (int j = 0; j < 4; j++) {
             out_vecs[i][j] = 0.0;
@@ -1012,13 +1011,13 @@ static void SetTexinfo_QuakeEd_New(
 
     //                        M[col][row]
     // S
-    out_vecs[0][axes.first] = M[0][0];
-    out_vecs[0][axes.second] = M[1][0];
+    out_vecs[0][axes.first] = M.at(0, 0);
+    out_vecs[0][axes.second] = M.at(0, 1);
     out_vecs[0][3] = shift[0];
 
     // T
-    out_vecs[1][axes.first] = M[0][1];
-    out_vecs[1][axes.second] = M[1][1];
+    out_vecs[1][axes.first] = M.at(1, 0);
+    out_vecs[1][axes.second] = M.at(1, 1);
     out_vecs[1][3] = shift[1];
 }
 
@@ -1783,7 +1782,7 @@ static void RotateMapFace(mapface_t *face, const qvec3d &angles)
     for (int i = 0; i < 2; i++) {
         const qvec4f in = texvecs.at(i);
         const qvec3f in_first3(in);
-
+        
         const qvec3f out_first3 = rotation * in_first3;
         newtexvecs[i] = {out_first3[0], out_first3[1], out_first3[2], in[3]};
     }
@@ -1810,6 +1809,7 @@ static void TranslateMapFace(mapface_t *face, const qvec3d &offset)
         // CHECK: precision loss here?
         out[3] += qv::dot(qvec3f(out), qvec3f(offset) * -1.0f);
         newtexvecs[i] = {out[0], out[1], out[2], out[3]};
+
     }
 
     face->set_texvecs(newtexvecs);
diff --git a/qbsp/portals.cc b/qbsp/portals.cc
index fc5fba5b..126773a9 100644
--- a/qbsp/portals.cc
+++ b/qbsp/portals.cc
@@ -382,7 +382,7 @@ static void MakeHeadnodePortals(const mapentity_t *entity, node_t *node)
                 pl->normal[i] = 1;
                 pl->dist = bounds[j][i];
             }
-            p->planenum = FindPlane(&pl->normal[0], pl->dist, &side);
+            p->planenum = FindPlane(pl->normal, pl->dist, &side);
 
             p->winding = BaseWindingForPlane(pl);
             if (side)
diff --git a/qbsp/qbsp.cc b/qbsp/qbsp.cc
index 4d17012e..94a8b752 100644
--- a/qbsp/qbsp.cc
+++ b/qbsp/qbsp.cc
@@ -130,7 +130,7 @@ static std::vector<std::tuple<size_t, face_t *>> AddBrushBevels(const brush_t *b
             auto flipped = map.planes[f->planenum];
             flipped.dist = -flipped.dist;
             VectorInverse(flipped.normal);
-            planenum = FindPlane(&flipped.normal[0], flipped.dist, nullptr);
+            planenum = FindPlane(flipped.normal, flipped.dist, nullptr);
         }
 
         int32_t outputplanenum = ExportMapPlane(planenum);
@@ -160,7 +160,7 @@ static std::vector<std::tuple<size_t, face_t *>> AddBrushBevels(const brush_t *b
                 else
                     new_plane.dist = -b->bounds.mins()[axis];
 
-                int32_t planenum = FindPlane(&new_plane.normal[0], new_plane.dist, nullptr);
+                int32_t planenum = FindPlane(new_plane.normal, new_plane.dist, nullptr);
                 int32_t outputplanenum = ExportMapPlane(planenum);
                 planes.emplace_back(outputplanenum, b->faces);
             }
@@ -248,7 +248,7 @@ static std::vector<std::tuple<size_t, face_t *>> AddBrushBevels(const brush_t *b
                         continue; // wasn't part of the outer hull
 
                     // add this plane
-                    int32_t planenum = FindPlane(&current.normal[0], current.dist, nullptr);
+                    int32_t planenum = FindPlane(current.normal, current.dist, nullptr);
                     int32_t outputplanenum = ExportMapPlane(planenum);
                     planes.emplace_back(outputplanenum, b->faces);
                 }
diff --git a/qbsp/test_qbsp.cc b/qbsp/test_qbsp.cc
index 1db6e6bb..87565ecf 100644
--- a/qbsp/test_qbsp.cc
+++ b/qbsp/test_qbsp.cc
@@ -164,7 +164,7 @@ TEST(qbsp, BrushMostlyOnSide1)
 {
     brush_t *brush = load128x128x32Brush();
 
-    vec3_t plane1normal = {-1, 0, 0};
+    qvec3d plane1normal {-1, 0, 0};
     vec_t plane1dist = -100;
 
     EXPECT_EQ(SIDE_FRONT, BrushMostlyOnSide(brush, plane1normal, plane1dist));
@@ -176,7 +176,7 @@ TEST(qbsp, BrushMostlyOnSide2)
 {
     brush_t *brush = load128x128x32Brush();
 
-    vec3_t plane1normal = {1, 0, 0};
+    qvec3d plane1normal {1, 0, 0};
     vec_t plane1dist = 100;
 
     EXPECT_EQ(SIDE_BACK, BrushMostlyOnSide(brush, plane1normal, plane1dist));
@@ -241,7 +241,7 @@ TEST(qbsp, SplitBrush)
 {
     brush_t *brush = load128x128x32Brush();
 
-    const vec3_t planenormal = {-1, 0, 0};
+    const qvec3d planenormal {-1, 0, 0};
     int planeside;
     const int planenum = FindPlane(planenormal, 0.0, &planeside);
 
@@ -282,7 +282,7 @@ TEST(qbsp, SplitBrushOnSide)
 {
     brush_t *brush = load128x128x32Brush();
 
-    const vec3_t planenormal = {-1, 0, 0};
+    const qvec3d planenormal {-1, 0, 0};
     int planeside;
     const int planenum = FindPlane(planenormal, -64.0, &planeside);
 
@@ -299,7 +299,7 @@ TEST(qbsp, SplitBrushOnSide)
 TEST(qbsp, MemLeaks) {
     brush_t *brush = load128x128x32Brush();
     
-    const vec3_t planenormal = { -1, 0, 0 };
+    const qvec3d planenormal { -1, 0, 0 };
     int planeside;
     const int planenum = FindPlane(planenormal, 0.0, &planeside);