diff --git a/core/src/geom.rs b/core/src/geom.rs
index d0b5208c..4c0a759d 100644
--- a/core/src/geom.rs
+++ b/core/src/geom.rs
@@ -1,6 +1,10 @@
 //! Basic geometric primitives.
 
-use crate::math::vec::{Vec2, Vec3};
+use crate::{
+    math::point::{Point2, Point3},
+    math::vec::{Vec2, Vec3},
+    render::Model,
+};
 
 pub use mesh::Mesh;
 
@@ -13,6 +17,12 @@ pub struct Vertex<P, A> {
     pub attrib: A,
 }
 
+/// Two-dimensional vertex type.
+pub type Vertex2<A, B = Model> = Vertex<Point2<B>, A>;
+
+/// Three-dimensional vertex type.
+pub type Vertex3<A, B = Model> = Vertex<Point3<B>, A>;
+
 /// Triangle, defined by three vertices.
 #[derive(Copy, Clone, Debug, Eq, PartialEq)]
 #[repr(transparent)]
diff --git a/core/src/geom/mesh.rs b/core/src/geom/mesh.rs
index ce204d5b..9082de1d 100644
--- a/core/src/geom/mesh.rs
+++ b/core/src/geom/mesh.rs
@@ -9,15 +9,12 @@ use alloc::{vec, vec::Vec};
 
 use crate::math::{
     mat::{Mat4x4, RealToReal},
+    point::Point3,
     space::Linear,
-    vec::Vec3,
 };
 use crate::render::Model;
 
-use super::{vertex, Normal3, Tri};
-
-/// Convenience type alias for a mesh vertex.
-pub type Vertex<A, B = Model> = super::Vertex<Vec3<B>, A>;
+use super::{vertex, Normal3, Tri, Vertex3};
 
 /// A triangle mesh.
 ///
@@ -32,7 +29,7 @@ pub struct Mesh<Attrib, Basis = Model> {
     /// to the `verts` vector. Several faces can share a vertex.
     pub faces: Vec<Tri<usize>>,
     /// The vertices of the mesh.
-    pub verts: Vec<Vertex<Attrib, Basis>>,
+    pub verts: Vec<Vertex3<Attrib, Basis>>,
 }
 
 /// A builder type for creating meshes.
@@ -54,16 +51,17 @@ impl<A, B> Mesh<A, B> {
     /// # Examples
     /// ```
     /// # use retrofire_core::geom::{Tri, Mesh, vertex};
-    /// # use retrofire_core::math::vec3;
+    /// # use retrofire_core::math::point::pt3;
     /// let verts = [
-    ///     vec3(0.0, 0.0, 0.0),
-    ///     vec3(1.0, 0.0, 0.0),
-    ///     vec3(0.0, 1.0, 0.0),
-    ///     vec3(0.0, 0.0, 1.0)
+    ///     pt3(0.0, 0.0, 0.0),
+    ///     pt3(1.0, 0.0, 0.0),
+    ///     pt3(0.0, 1.0, 0.0),
+    ///     pt3(0.0, 0.0, 1.0)
     /// ]
     /// .map(|v| vertex(v, ()));
     ///
     /// let faces = [
+    ///     // Indices point to the verts array
     ///     Tri([0, 1, 2]),
     ///     Tri([0, 1, 3]),
     ///     Tri([0, 2, 3]),
@@ -78,7 +76,7 @@ impl<A, B> Mesh<A, B> {
     pub fn new<F, V>(faces: F, verts: V) -> Self
     where
         F: IntoIterator<Item = Tri<usize>>,
-        V: IntoIterator<Item = Vertex<A, B>>,
+        V: IntoIterator<Item = Vertex3<A, B>>,
     {
         let faces: Vec<_> = faces.into_iter().collect();
         let verts: Vec<_> = verts.into_iter().collect();
@@ -129,16 +127,16 @@ impl<A> Builder<A> {
     }
 
     /// Appends a vertex with the given position and attribute.
-    pub fn push_vert(&mut self, pos: Vec3, attrib: A) {
+    pub fn push_vert(&mut self, pos: Point3, attrib: A) {
         self.mesh.verts.push(vertex(pos.to(), attrib));
     }
 
     /// Appends all the vertices yielded by the given iterator.
     pub fn push_verts<Vs>(&mut self, verts: Vs)
     where
-        Vs: IntoIterator<Item = (Vec3, A)>,
+        Vs: IntoIterator<Item = (Point3, A)>,
     {
-        let vs = verts.into_iter().map(|(v, a)| vertex(v.to(), a));
+        let vs = verts.into_iter().map(|(p, a)| vertex(p.to(), a));
         self.mesh.verts.extend(vs);
     }
 
@@ -167,7 +165,7 @@ impl Builder<()> {
                 .mesh
                 .verts
                 .into_iter()
-                .map(|v| vertex(tf.apply(&v.pos), v.attrib))
+                .map(|v| vertex(tf.apply_pt(&v.pos), v.attrib))
                 .collect(),
         };
         mesh.into_builder()
@@ -257,6 +255,7 @@ mod tests {
     use core::f32::consts::FRAC_1_SQRT_2;
 
     use crate::geom::vertex;
+    use crate::math::point::pt3;
     use crate::math::vec3;
     use crate::prelude::splat;
 
@@ -268,9 +267,9 @@ mod tests {
         let _: Mesh<()> = Mesh::new(
             [Tri([0, 1, 2]), Tri([1, 2, 3])],
             [
-                vertex(vec3(0.0, 0.0, 0.0), ()),
-                vertex(vec3(1.0, 1.0, 1.0), ()),
-                vertex(vec3(2.0, 2.0, 2.0), ()),
+                vertex(pt3(0.0, 0.0, 0.0), ()),
+                vertex(pt3(1.0, 1.0, 1.0), ()),
+                vertex(pt3(2.0, 2.0, 2.0), ()),
             ],
         );
     }
@@ -281,9 +280,9 @@ mod tests {
         let mut b = Mesh::builder();
         b.push_faces([[0, 1, 2], [1, 2, 3]]);
         b.push_verts([
-            (vec3(0.0, 0.0, 0.0), ()),
-            (vec3(1.0, 1.0, 1.0), ()),
-            (vec3(2.0, 2.0, 2.0), ()),
+            (pt3(0.0, 0.0, 0.0), ()),
+            (pt3(1.0, 1.0, 1.0), ()),
+            (pt3(2.0, 2.0, 2.0), ()),
         ]);
         _ = b.build();
     }
@@ -295,10 +294,10 @@ mod tests {
         let mut b = Mesh::builder();
         b.push_faces([[0, 2, 1], [0, 1, 3], [0, 3, 2]]);
         b.push_verts([
-            (vec3(0.0, 0.0, 0.0), ()),
-            (vec3(1.0, 0.0, 0.0), ()),
-            (vec3(0.0, 1.0, 0.0), ()),
-            (vec3(0.0, 0.0, 1.0), ()),
+            (pt3(0.0, 0.0, 0.0), ()),
+            (pt3(1.0, 0.0, 0.0), ()),
+            (pt3(0.0, 1.0, 0.0), ()),
+            (pt3(0.0, 0.0, 1.0), ()),
         ]);
         let b = b.with_vertex_normals();
 
diff --git a/core/src/math.rs b/core/src/math.rs
index c3e41915..557969ee 100644
--- a/core/src/math.rs
+++ b/core/src/math.rs
@@ -30,8 +30,22 @@ pub mod approx;
 pub mod color;
 pub mod float;
 pub mod mat;
+pub mod point;
 pub mod rand;
 pub mod space;
 pub mod spline;
 pub mod vary;
 pub mod vec;
+
+pub trait Lerp {
+    fn lerp(&self, other: &Self, t: f32) -> Self;
+}
+
+impl<T> Lerp for T
+where
+    T: Affine<Diff: Linear<Scalar = f32>>,
+{
+    fn lerp(&self, other: &Self, t: f32) -> Self {
+        self.add(&other.sub(self).mul(t))
+    }
+}
diff --git a/core/src/math/angle.rs b/core/src/math/angle.rs
index b823a7bb..5a8a330e 100644
--- a/core/src/math/angle.rs
+++ b/core/src/math/angle.rs
@@ -10,6 +10,7 @@ use crate::math::vec::Vector;
 
 #[cfg(feature = "fp")]
 use crate::math::float::f32;
+use crate::math::vary::ZDiv;
 #[cfg(feature = "fp")]
 use crate::math::vec::{vec2, vec3, Vec2, Vec3};
 
@@ -463,6 +464,8 @@ impl Linear for Angle {
     }
 }
 
+impl ZDiv for Angle {}
+
 //
 // Foreign trait impls
 //
@@ -571,7 +574,6 @@ mod tests {
     use core::f32::consts::{PI, TAU};
 
     use crate::assert_approx_eq;
-    use crate::math::vary::Vary;
 
     use super::*;
 
@@ -678,6 +680,7 @@ mod tests {
 
     #[test]
     fn varying() {
+        use crate::math::vary::Vary;
         let mut i = degs(45.0).vary(degs(15.0), Some(4));
 
         assert_approx_eq!(i.next(), Some(degs(45.0)));
diff --git a/core/src/math/color.rs b/core/src/math/color.rs
index 3c8c5850..59af4840 100644
--- a/core/src/math/color.rs
+++ b/core/src/math/color.rs
@@ -1,13 +1,18 @@
 //! Colors and color spaces.
 
-use core::array;
-use core::fmt::{self, Debug, Formatter};
-use core::marker::PhantomData;
-use core::ops::Index;
-
-use crate::math::float::f32;
-use crate::math::space::{Affine, Linear};
-use crate::math::vec::Vector;
+use core::{
+    array,
+    fmt::{self, Debug, Formatter},
+    marker::PhantomData,
+    ops::Index,
+};
+
+use crate::math::{
+    float::f32,
+    space::{Affine, Linear},
+    vary::ZDiv,
+    vec::Vector,
+};
 
 //
 // Types
@@ -531,6 +536,8 @@ impl<Sp, const DIM: usize> Linear for Color<[f32; DIM], Sp> {
     }
 }
 
+impl<Sc, Sp, const N: usize> ZDiv for Color<[Sc; N], Sp> where Sc: ZDiv + Copy {}
+
 //
 // Foreign trait impls
 //
diff --git a/core/src/math/mat.rs b/core/src/math/mat.rs
index 9a19e2e4..1a86a8bc 100644
--- a/core/src/math/mat.rs
+++ b/core/src/math/mat.rs
@@ -1,16 +1,23 @@
 #![allow(clippy::needless_range_loop)]
 
 //! Matrices and linear transforms.
+//!
+//! TODO Docs
 
-use core::array;
-use core::fmt::{self, Debug, Formatter};
-use core::marker::PhantomData;
-use core::ops::Range;
+use core::{
+    array,
+    fmt::{self, Debug, Formatter},
+    marker::PhantomData,
+    ops::Range,
+};
 
 use crate::render::{NdcToScreen, ViewToProj};
 
-use super::space::{Linear, Proj4, Real};
-use super::vec::{ProjVec4, Vec2, Vec2u, Vec3, Vector};
+use super::{
+    point::{Point2, Point3},
+    space::{Linear, Proj4, Real},
+    vec::{ProjVec4, Vec2, Vec2u, Vec3, Vector},
+};
 
 /// A linear transform from one space (or basis) to another.
 ///
@@ -199,9 +206,16 @@ impl<Src, Dst> Mat3x3<RealToReal<2, Src, Dst>> {
     /// ```
     #[must_use]
     pub fn apply(&self, v: &Vec2<Src>) -> Vec2<Dst> {
-        let v = [v.x(), v.y(), 1.0].into();
+        let v = [v.x(), v.y(), 1.0].into(); // TODO w=0.0
         array::from_fn(|i| self.row_vec(i).dot(&v)).into()
     }
+
+    // TODO Add trait to overload apply or similar
+    #[must_use]
+    pub fn apply_pt(&self, p: &Point2<Src>) -> Point2<Dst> {
+        let p = [p.x(), p.y(), 1.0].into();
+        array::from_fn(|i| self.row_vec(i).dot(&p)).into()
+    }
 }
 
 impl<Src, Dst> Mat4x4<RealToReal<3, Src, Dst>> {
@@ -218,10 +232,17 @@ impl<Src, Dst> Mat4x4<RealToReal<3, Src, Dst>> {
     /// ```
     #[must_use]
     pub fn apply(&self, v: &Vec3<Src>) -> Vec3<Dst> {
-        let v = [v.x(), v.y(), v.z(), 1.0].into();
+        let v = [v.x(), v.y(), v.z(), 1.0].into(); // TODO w=0.0
         array::from_fn(|i| self.row_vec(i).dot(&v)).into()
     }
 
+    // TODO Add trait to overload apply or similar
+    #[must_use]
+    pub fn apply_pt(&self, p: &Point3<Src>) -> Point3<Dst> {
+        let p = [p.x(), p.y(), p.z(), 1.0].into();
+        array::from_fn(|i| self.row_vec(i).dot(&p)).into()
+    }
+
     /// Returns the determinant of `self`.
     ///
     /// Given a matrix M,
@@ -365,8 +386,8 @@ impl<Src> Mat4x4<RealToProj<Src>> {
     ///         \ ·  ·  M33 / \  1 /
     /// ```
     #[must_use]
-    pub fn apply(&self, v: &Vec3<Src>) -> ProjVec4 {
-        let v = Vector::from([v.x(), v.y(), v.z(), 1.0]);
+    pub fn apply(&self, p: &Point3<Src>) -> ProjVec4 {
+        let v = Vector::from([p.x(), p.y(), p.z(), 1.0]);
         [
             self.row_vec(0).dot(&v),
             self.row_vec(1).dot(&v),
@@ -602,13 +623,17 @@ pub fn perspective(
 /// # Parameters
 /// * `lbn`: The left-bottom-near corner of the projection box.
 /// * `rtf`: The right-bottom-far corner of the projection box.
-pub fn orthographic(lbn: Vec3, rtf: Vec3) -> Mat4x4<ViewToProj> {
-    let [dx, dy, dz] = (rtf - lbn).0;
-    let [sx, sy, sz] = (rtf + lbn).0;
+// TODO Change to take points
+pub fn orthographic(lbn: Point3, rtf: Point3) -> Mat4x4<ViewToProj> {
+    let half_d = 0.5 * (rtf - lbn);
+    let center_pt = lbn + half_d;
+
+    let [dx, dy, dz] = half_d.0;
+    let [cx, cy, cz] = center_pt.0;
     [
-        [2.0 / dx, 0.0, 0.0, -sx / dx],
-        [0.0, 2.0 / dy, 0.0, -sy / dy],
-        [0.0, 0.0, 2.0 / dz, -sz / dz],
+        [1.0 / dx, 0.0, 0.0, -cx / dx],
+        [0.0, 1.0 / dy, 0.0, -cy / dy],
+        [0.0, 0.0, 1.0 / dz, -cz / dz],
         [0.0, 0.0, 0.0, 1.0],
     ]
     .into()
@@ -639,6 +664,7 @@ mod tests {
 
     #[cfg(feature = "fp")]
     use crate::math::angle::degs;
+    use crate::math::point::pt3;
 
     use super::*;
 
@@ -652,6 +678,7 @@ mod tests {
 
     mod mat3x3 {
         use super::*;
+        use crate::math::point::pt2;
 
         const MAT: Mat3x3<Map> = Matrix::new([
             [0.0, 1.0, 2.0], //
@@ -696,6 +723,7 @@ mod tests {
                 [0.0, 0.0, 1.0],
             ]);
             assert_eq!(m.apply(&vec2(1.0, 2.0)), vec2(2.0, -6.0));
+            assert_eq!(m.apply_pt(&pt2(2.0, -1.0)), pt2(4.0, 3.0));
         }
 
         #[test]
@@ -706,6 +734,7 @@ mod tests {
                 [0.0, 0.0, 1.0],
             ]);
             assert_eq!(m.apply(&vec2(1.0, 2.0)), vec2(3.0, -1.0));
+            assert_eq!(m.apply_pt(&pt2(2.0, -1.0)), pt2(4.0, -4.0));
         }
 
         #[test]
@@ -723,6 +752,7 @@ mod tests {
 
     mod mat4x4 {
         use super::*;
+        use crate::math::point::pt3;
 
         const MAT: Mat4x4<Map> = Matrix::new([
             [0.0, 1.0, 2.0, 3.0],
@@ -753,17 +783,25 @@ mod tests {
         }
 
         #[test]
-        fn scaling_vec3() {
+        fn scaling() {
             let m = scale(vec3(1.0, -2.0, 3.0));
+
             let v = vec3(0.0, 4.0, -3.0);
             assert_eq!(m.apply(&v), vec3(0.0, -8.0, -9.0));
+
+            let p = pt3(4.0, 0.0, -3.0);
+            assert_eq!(m.apply_pt(&p), pt3(4.0, 0.0, -9.0));
         }
 
         #[test]
-        fn translation_vec3() {
+        fn translation() {
             let m = translate(vec3(1.0, 2.0, 3.0));
+
             let v = vec3(0.0, 5.0, -3.0);
             assert_eq!(m.apply(&v), vec3(1.0, 7.0, 0.0));
+
+            let p = pt3(3.0, 5.0, 0.0);
+            assert_eq!(m.apply_pt(&p), pt3(4.0, 7.0, 3.0));
         }
 
         #[cfg(feature = "fp")]
@@ -775,6 +813,10 @@ mod tests {
                 m.apply(&vec3(0.0, 0.0, 1.0)),
                 vec3(0.0, 1.0, 0.0)
             );
+            assert_approx_eq!(
+                m.apply_pt(&pt3(0.0, -2.0, 0.0)),
+                pt3(0.0, 0.0, 2.0)
+            );
         }
 
         #[cfg(feature = "fp")]
@@ -786,6 +828,10 @@ mod tests {
                 m.apply(&vec3(1.0, 0.0, 0.0)),
                 vec3(0.0, 0.0, 1.0)
             );
+            assert_approx_eq!(
+                m.apply_pt(&pt3(0.0, 0.0, -2.0)),
+                pt3(2.0, 0.0, 0.0)
+            );
         }
 
         #[cfg(feature = "fp")]
@@ -797,6 +843,10 @@ mod tests {
                 m.apply(&vec3(0.0, 1.0, 0.0)),
                 vec3(1.0, 0.0, 0.0)
             );
+            assert_approx_eq!(
+                m.apply_pt(&pt3(-2.0, 0.0, 0.0)),
+                pt3(0.0, 2.0, 0.0)
+            );
         }
 
         #[test]
@@ -884,19 +934,19 @@ mod tests {
 
     #[test]
     fn orthographic_box_maps_to_unit_cube() {
-        let lbn = vec3(-20.0, 0.0, 0.01);
-        let rtf = vec3(100.0, 50.0, 100.0);
+        let lbn = pt3(-20.0, 0.0, 0.01);
+        let rtf = pt3(100.0, 50.0, 100.0);
 
         let m = orthographic(lbn, rtf);
 
         assert_approx_eq!(m.apply(&lbn.to()), [-1.0, -1.0, -1.0, 1.0].into());
-        assert_approx_eq!(m.apply(&rtf.to()), splat(1.0));
+        assert_approx_eq!(m.apply(&rtf.to()), [1.0, 1.0, 1.0, 1.0].into());
     }
 
     #[test]
     fn perspective_frustum_maps_to_unit_cube() {
-        let left_bot_near = vec3(-0.125, -0.0625, 0.1);
-        let right_top_far = vec3(125.0, 62.5, 100.0);
+        let left_bot_near = pt3(-0.125, -0.0625, 0.1);
+        let right_top_far = pt3(125.0, 62.5, 100.0);
 
         let m = perspective(0.8, 2.0, 0.1..100.0);
 
@@ -904,6 +954,6 @@ mod tests {
         assert_approx_eq!(lbn / lbn.w(), [-1.0, -1.0, -1.0, 1.0].into());
 
         let rtf = m.apply(&right_top_far);
-        assert_approx_eq!(rtf / rtf.w(), splat(1.0));
+        assert_approx_eq!(rtf / rtf.w(), [1.0, 1.0, 1.0, 1.0].into());
     }
 }
diff --git a/core/src/math/point.rs b/core/src/math/point.rs
new file mode 100644
index 00000000..c3e0de55
--- /dev/null
+++ b/core/src/math/point.rs
@@ -0,0 +1,236 @@
+use core::{
+    array,
+    fmt::{Debug, Formatter},
+    marker::PhantomData as Pd,
+    ops::{Add, Index, Sub},
+};
+
+use crate::math::vary::ZDiv;
+use crate::math::{
+    space::{Affine, Linear, Real},
+    vec::Vector,
+    ApproxEq,
+};
+
+#[repr(transparent)]
+pub struct Point<Repr, Space = ()>(pub Repr, Pd<Space>);
+
+/// A 2-point with `f32` components.
+pub type Point2<Basis = ()> = Point<[f32; 2], Real<2, Basis>>;
+/// A 3-point with `f32` components.
+pub type Point3<Basis = ()> = Point<[f32; 3], Real<3, Basis>>;
+
+/// A 2-point with `u32` components.
+pub type Point2u<Basis = ()> = Point<[u32; 2], Real<2, Basis>>;
+
+/// Returns a real 2-point with `x` and `y` components.
+pub const fn pt2<Sc, B>(x: Sc, y: Sc) -> Point<[Sc; 2], Real<2, B>> {
+    Point([x, y], Pd)
+}
+/// Returns a real 3-point with `x`, `y`, and `z` components.
+pub const fn pt3<Sc, B>(x: Sc, y: Sc, z: Sc) -> Point<[Sc; 3], Real<3, B>> {
+    Point([x, y, z], Pd)
+}
+
+impl<R, Sp> Point<R, Sp> {
+    #[inline]
+    pub fn new(repr: R) -> Self {
+        Self(repr, Pd)
+    }
+
+    /// Returns a point with value equal to `self` but in space `S`.
+    // TODO Cannot be const (yet?) due to E0493 :(
+    #[inline]
+    pub fn to<S>(self) -> Point<R, S> {
+        Point(self.0, Pd)
+    }
+
+    /// Returns the vector equivalent to `self`.
+    // TODO Cannot be const (yet?) due to E0493 :(
+    #[inline]
+    pub fn to_vec(self) -> Vector<R, Sp> {
+        Vector::new(self.0)
+    }
+}
+
+impl<Sc, Sp, const N: usize> Point<[Sc; N], Sp> {
+    /// Returns a vector of the same dimension as `self` by applying `f`
+    /// component-wise.
+    #[inline]
+    #[must_use]
+    pub fn map<T>(self, f: impl FnMut(Sc) -> T) -> Point<[T; N], Sp> {
+        self.0.map(f).into()
+    }
+}
+
+impl<const N: usize, B> Point<[f32; N], Real<N, B>> {
+    #[cfg(feature = "fp")]
+    #[inline]
+    pub fn distance(&self, other: &Self) -> f32 {
+        Affine::sub(self, other).len()
+    }
+    #[inline]
+    pub fn distance_sqr(&self, other: &Self) -> f32 {
+        Affine::sub(self, other).len_sqr()
+    }
+
+    /// Returns `self` clamped component-wise to the given range.
+    ///
+    /// In other words, for each component `self[i]`, the result `r` has
+    /// `r[i]` equal to `self[i].clamp(min[i], max[i])`.
+    ///
+    /// # Examples
+    /// ```
+    /// # use retrofire_core::math::point::{pt3, Point3};
+    /// let pt = pt3::<f32, ()>(0.5, 1.5, -2.0);
+    /// // Clamp to the unit cube
+    /// let clamped = pt.clamp(&pt3(0.0, 0.0, 0.0), &pt3(1.0, 1.0, 1.0));
+    /// assert_eq!(clamped, pt3(0.5, 1.0, -0.0));
+    #[must_use]
+    pub fn clamp(&self, min: &Self, max: &Self) -> Self {
+        array::from_fn(|i| self.0[i].clamp(min.0[i], max.0[i])).into()
+    }
+}
+
+impl<R, B, Sc> Point<R, Real<2, B>>
+where
+    R: Index<usize, Output = Sc>,
+    Sc: Copy,
+{
+    /// Returns the x component of `self`.
+    #[inline]
+    pub fn x(&self) -> Sc {
+        self.0[0]
+    }
+    /// Returns the y component of `self`.
+    #[inline]
+    pub fn y(&self) -> Sc {
+        self.0[1]
+    }
+}
+
+impl<R, Sc, B> Point<R, Real<3, B>>
+where
+    R: Index<usize, Output = Sc>,
+    Sc: Copy,
+{
+    /// Returns the x component of `self`.
+    #[inline]
+    pub fn x(&self) -> Sc {
+        self.0[0]
+    }
+    /// Returns the y component of `self`.
+    #[inline]
+    pub fn y(&self) -> Sc {
+        self.0[1]
+    }
+    /// Returns the z component of `self`.
+    #[inline]
+    pub fn z(&self) -> Sc {
+        self.0[2]
+    }
+}
+
+impl<Sc, Sp, const N: usize> Affine for Point<[Sc; N], Sp>
+where
+    Sc: Linear<Scalar = Sc> + Copy,
+{
+    type Space = Sp;
+    type Diff = Vector<[Sc; N], Sp>;
+    const DIM: usize = N;
+
+    #[inline]
+    fn add(&self, other: &Self::Diff) -> Self {
+        // TODO Profile performance of array::from_fn
+        Self(array::from_fn(|i| self.0[i].add(&other.0[i])), Pd)
+    }
+    #[inline]
+    fn sub(&self, other: &Self) -> Self::Diff {
+        Vector::new(array::from_fn(|i| self.0[i].sub(&other.0[i])))
+    }
+}
+
+impl<Sc, Sp, const N: usize> ZDiv for Point<[Sc; N], Sp>
+where
+    Sc: ZDiv + Copy,
+{
+    fn z_div(self, z: f32) -> Self {
+        Self(self.0.map(|c| c.z_div(z)), Pd)
+    }
+}
+
+impl<Sc: ApproxEq, Sp, const N: usize> ApproxEq<Self, Sc>
+    for Point<[Sc; N], Sp>
+{
+    fn approx_eq_eps(&self, other: &Self, eps: &Sc) -> bool {
+        self.0.approx_eq_eps(&other.0, eps)
+    }
+    fn relative_epsilon() -> Sc {
+        Sc::relative_epsilon()
+    }
+}
+
+//
+// Foreign trait impls
+//
+
+// Manual impls of Copy, Clone, Eq, and PartialEq to avoid
+// superfluous where S: Trait bound
+
+impl<R: Copy, S> Copy for Point<R, S> {}
+
+impl<R: Clone, S> Clone for Point<R, S> {
+    fn clone(&self) -> Self {
+        Self(self.0.clone(), Pd)
+    }
+}
+
+impl<R: Default, S> Default for Point<R, S> {
+    fn default() -> Self {
+        Self(R::default(), Pd)
+    }
+}
+
+impl<R: Debug, Sp: Debug + Default> Debug for Point<R, Sp> {
+    fn fmt(&self, f: &mut Formatter<'_>) -> core::fmt::Result {
+        write!(f, "Point<{:?}>", Sp::default())?;
+        Debug::fmt(&self.0, f)
+    }
+}
+
+impl<R: Eq, S> Eq for Point<R, S> {}
+
+impl<R: PartialEq, S> PartialEq for Point<R, S> {
+    fn eq(&self, other: &Self) -> bool {
+        self.0 == other.0
+    }
+}
+
+impl<R, Sp> From<R> for Point<R, Sp> {
+    #[inline]
+    fn from(repr: R) -> Self {
+        Self(repr, Pd)
+    }
+}
+
+impl<R, Sp> Add<<Self as Affine>::Diff> for Point<R, Sp>
+where
+    Self: Affine,
+{
+    type Output = Self;
+
+    fn add(self, other: <Self as Affine>::Diff) -> Self {
+        Affine::add(&self, &other)
+    }
+}
+
+impl<R, Sp> Sub for Point<R, Sp>
+where
+    Self: Affine,
+{
+    type Output = <Self as Affine>::Diff;
+
+    fn sub(self, other: Self) -> Self::Output {
+        Affine::sub(&self, &other)
+    }
+}
diff --git a/core/src/math/space.rs b/core/src/math/space.rs
index 7e37f0d8..e141e8a1 100644
--- a/core/src/math/space.rs
+++ b/core/src/math/space.rs
@@ -2,9 +2,10 @@
 //!
 //! TODO
 
+use core::fmt::{Debug, Formatter};
 use core::marker::PhantomData;
 
-use crate::math::vary::{Iter, Vary};
+use crate::math::vary::{Iter, Vary, ZDiv};
 
 /// Trait for types representing elements of an affine space.
 ///
@@ -29,6 +30,25 @@ pub trait Affine: Sized {
     ///
     /// `sub` is anti-commutative: `v.sub(w) == w.sub(v).neg()`.
     fn sub(&self, other: &Self) -> Self::Diff;
+
+    /// Linearly interpolates between `self` and `other`.
+    ///
+    /// This method does not panic if `t < 0.0` or `t > 1.0`, or if `t`
+    /// is a `NaN`, but the return value in those cases is unspecified.
+    /// Individual implementations may offer stronger guarantees.
+    ///
+    /// # Examples
+    /// ```
+    /// # use retrofire_core::math::space::Affine;
+    /// assert_eq!(2.0.lerp(&5.0, 0.0), 2.0);
+    /// assert_eq!(2.0.lerp(&5.0, 0.5), 3.5);
+    /// assert_eq!(2.0.lerp(&5.0, 1.0), 5.0);
+    /// ```
+    // TODO Now duplicated in Vary
+    #[inline]
+    fn lerp(&self, other: &Self, t: <Self::Diff as Linear>::Scalar) -> Self {
+        self.add(&other.sub(self).mul(t))
+    }
 }
 
 /// Trait for types representing elements of a linear space (vector space).
@@ -67,12 +87,14 @@ pub trait Linear: Affine<Diff = Self> {
     fn mul(&self, scalar: Self::Scalar) -> Self;
 }
 
-/// Tag type for real vector spaces (Euclidean spaces) of dimension `DIM`.
+/// Tag type for real vector spaces and Euclidean spaces.
+///
 /// For example, the type `Real<3>` corresponds to ℝ³.
 #[derive(Copy, Clone, Default, Eq, PartialEq)]
 pub struct Real<const DIM: usize, Basis = ()>(PhantomData<Basis>);
 
 /// Tag type for the projective 4-space over reals, 𝗣<sub>4</sub>(ℝ).
+///
 /// The properties of this space make it useful for implementing perspective
 /// projection. Clipping is also done in the projective space.
 #[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
@@ -155,7 +177,7 @@ impl Affine for u32 {
 
 impl<V: Clone> Vary for V
 where
-    Self: Linear<Scalar = f32>,
+    Self: Affine<Diff: Linear<Scalar = f32> + Clone> + ZDiv,
 {
     type Iter = Iter<Self>;
     type Diff = <Self as Affine>::Diff;
@@ -174,9 +196,20 @@ where
     fn step(&self, delta: &Self::Diff) -> Self {
         self.add(delta)
     }
+}
 
-    fn z_div(&self, z: f32) -> Self {
-        self.mul(z.recip())
+impl<const DIM: usize, B> Debug for Real<DIM, B>
+where
+    B: Debug + Default,
+{
+    fn fmt(&self, f: &mut Formatter<'_>) -> core::fmt::Result {
+        const DIMS: [&str; 4] = ["", "²", "³", "⁴"];
+        let b = B::default();
+        if let Some(dim) = DIMS.get(DIM - 1) {
+            write!(f, "ℝ{dim}<{b:?}>")
+        } else {
+            write!(f, "ℝ^{DIM}<{b:?}>")
+        }
     }
 }
 
diff --git a/core/src/math/spline.rs b/core/src/math/spline.rs
index 796e1755..35867b1b 100644
--- a/core/src/math/spline.rs
+++ b/core/src/math/spline.rs
@@ -1,11 +1,13 @@
-//! Bezier curves and splines.
+//! Bézier curves and splines.
 
 use alloc::vec::Vec;
 
-use crate::math::space::Linear;
-use crate::math::Vary;
+use crate::math::{
+    space::{Affine, Linear},
+    vary::ZDiv,
+};
 
-/// A cubic Bezier curve, defined by four control points.
+/// A cubic Bézier curve, defined by four control points.
 ///
 /// TODO More info about Beziers
 ///
@@ -56,9 +58,7 @@ where
 
 impl<T> CubicBezier<T>
 where
-    // T: Affine + Clone,
-    // T::Diff: Linear<Scalar = f32> + Clone,
-    T: Linear<Scalar = f32> + Clone,
+    T: Affine<Diff: Linear<Scalar = f32>> + Clone,
 {
     /// Evaluates the value of `self` at `t`
     ///
@@ -75,54 +75,74 @@ where
 }
 impl<T> CubicBezier<T>
 where
-    T: Linear<Scalar = f32> + Clone,
+    T: Affine<Diff: Linear<Scalar = f32>> + Clone,
 {
     pub fn fast_eval(&self, t: f32) -> T {
         let [p0, p1, p2, p3] = &self.0;
         step(t, p0, p3, |t| {
+            // Rewrite the parametric equation into a form where three of the
+            // coefficients are vectors, their linear combination added to `p0`
+            // so the equation can be expressed for affine types:
+            //
             //   (p3 - p0) * t^3 + (p1 - p2) * 3t^3
             // + (p0 + p2) * 3t^2 - p1 * 6t^2
             // + (p1 - p0) * 3t
             // + p0
-
-            let term3 = &p1.sub(p2).mul(3.0).add(&p3.sub(p0)).mul(t);
-            let term2 = &p1.mul(-2.0).add(p0).add(p2).mul(3.0);
-            let term1 = &p1.sub(p0).mul(3.0 * t);
-            let term0 = p0;
-
-            term3.add(term2).mul(t * t).add(term1).add(term0)
+            // = ((p3 - p0 + 3(p1 - p2)) * t^3
+            // + 3(p0 - p1 + p2 - p1) * t^2
+            // + 3(p1 - p0) * t
+            // + p0
+            // = ((((p3 - p0 + 3(p1 - p2))) * t
+            //      + 3(p0 - p1 + p2 - p1)) * t)
+            //          + 3(p1 - p0)) * t)
+            //              + p0
+
+            let p1_p0_3: T::Diff = p1.sub(p0).mul(3.0);
+            let p1_p2_3: T::Diff = p1.sub(p2).mul(3.0);
+
+            let co3: T::Diff = p3.sub(p0).add(&p1_p2_3);
+            let co2: T::Diff = p1_p0_3.add(&p1_p2_3);
+            let co1: T::Diff = p1_p0_3;
+            let co0: &T = p0;
+
+            // Add a linear combination of the three vectors to `p0`
+            // to get the result:
+            co0.add(&co3.mul(t).sub(&co2).mul(t).add(&co1).mul(t))
         })
     }
 
     /// Returns the tangent, or direction vector, of `self` at `t`.
     ///
     /// Clamps `t` to the range [0, 1].
-    pub fn tangent(&self, t: f32) -> T {
+    pub fn tangent(&self, t: f32) -> T::Diff {
         let [p0, p1, p2, p3] = &self.0;
         let t = t.clamp(0.0, 1.0);
 
-        //   (p1 - p2) * 9t^2 + (p3 - p0) * 3t^2
-        // + (p0 + p2) * 6t - p1 * 12t
-        // + (p1 - p0) * 3
+        //   3 (3 (p1 - p2) + (p3 - p0)) * t^2
+        // + 6 ((p0 - p1 + p2 - p1) * t
+        // + 3 (p1 - p0)
 
-        let term2 = p1.sub(p2).mul(3.0).add(&p3.sub(p0)).mul(t * t);
-        let term1 = p1.mul(-2.0).add(p0).add(p2).mul(2.0 * t);
-        let term0 = p1.sub(p0);
+        let co2: T::Diff = p1.sub(p2).mul(3.0).add(&p3.sub(p0));
+        let co1: T::Diff = p0.sub(p1).add(&p2.sub(p1)).mul(2.0);
+        let co0: T::Diff = p1.sub(p0);
 
-        term2.add(&term1).add(&term0).mul(3.0)
+        co2.mul(t).add(&co1).mul(t).add(&co0).mul(3.0)
     }
 }
 
 /// A curve composed of one or more concatenated
-/// [cubic Bezier curves][CubicBezier].
+/// [cubic Bézier curves][CubicBezier].
 #[derive(Debug, Clone, Eq, PartialEq)]
 pub struct BezierSpline<T>(Vec<T>);
 
-impl<T: Linear<Scalar = f32> + Copy> BezierSpline<T> {
-    /// Creates a bezier curve from the given control points. The number of
+impl<T> BezierSpline<T>
+where
+    T: Affine<Diff: Linear<Scalar = f32> + Copy> + ZDiv + Copy,
+{
+    /// Creates a Bézier curve from the given control points. The number of
     /// elements in `pts` must be 3n + 1 for some positive integer n.
     ///
-    /// Consecutive points in `pts` make up Bezier curves such that:
+    /// Consecutive points in `pts` make up Bézier curves such that:
     /// * `pts[0..=3]` define the first curve,
     /// * `pts[3..=6]` define the second curve,
     ///
@@ -153,7 +173,7 @@ impl<T: Linear<Scalar = f32> + Copy> BezierSpline<T> {
     /// Returns the tangent of `self` at `t`.
     ///
     /// Clamps `t` to the range [0, 1].
-    pub fn tangent(&self, t: f32) -> T {
+    pub fn tangent(&self, t: f32) -> T::Diff {
         let (t, seg) = self.segment(t);
         CubicBezier(seg).tangent(t)
     }
@@ -201,7 +221,7 @@ impl<T: Linear<Scalar = f32> + Copy> BezierSpline<T> {
     /// let approx = curve.approximate(|err| err.len_sqr() < 0.01*0.01);
     /// assert_eq!(approx.len(), 17);
     /// ```
-    pub fn approximate(&self, halt: impl Fn(&T) -> bool) -> Vec<T> {
+    pub fn approximate(&self, halt: impl Fn(&T::Diff) -> bool) -> Vec<T> {
         let len = self.0.len();
         let mut res = Vec::with_capacity(3 * len);
         self.do_approx(0.0, 1.0, 10 + len.ilog2(), &halt, &mut res);
@@ -214,7 +234,7 @@ impl<T: Linear<Scalar = f32> + Copy> BezierSpline<T> {
         a: f32,
         b: f32,
         max_dep: u32,
-        halt: &impl Fn(&T) -> bool,
+        halt: &impl Fn(&T::Diff) -> bool,
         accum: &mut Vec<T>,
     ) {
         let mid = a.lerp(&b, 0.5);
@@ -270,7 +290,7 @@ mod tests {
     }
 
     #[test]
-    fn bezier_spline_eval_eq_fasteval() {
+    fn bezier_spline_eval_eq_fast_eval() {
         let b: CubicBezier<Vec2> = CubicBezier(
             [[0.0, 0.0], [0.0, 2.0], [1.0, -1.0], [1.0, 1.0]].map(Vec2::from),
         );
diff --git a/core/src/math/vary.rs b/core/src/math/vary.rs
index 5474d4ce..4f37dab0 100644
--- a/core/src/math/vary.rs
+++ b/core/src/math/vary.rs
@@ -5,13 +5,20 @@
 
 use core::mem;
 
+pub trait ZDiv: Sized {
+    #[must_use]
+    fn z_div(self, _z: f32) -> Self {
+        self
+    }
+}
+
 /// A trait for types that can be linearly interpolated and distributed
 /// between two endpoints.
 ///
 /// This trait is designed particularly for *varyings:* types that are
 /// meant to be interpolated across the face of a polygon when rendering,
 /// but the methods are useful for various purposes.
-pub trait Vary: Sized + Clone {
+pub trait Vary: ZDiv + Sized + Clone {
     /// The iterator returned by the [vary][Self::vary] method.
     type Iter: Iterator<Item = Self>;
     /// The difference type of `Self`.
@@ -50,10 +57,6 @@ pub trait Vary: Sized + Clone {
     #[must_use]
     fn step(&self, delta: &Self::Diff) -> Self;
 
-    /// Performs perspective division.
-    #[must_use]
-    fn z_div(&self, z: f32) -> Self;
-
     /// Linearly interpolates between `self` and `other`.
     ///
     /// This method does not panic if `t < 0.0` or `t > 1.0`, or if `t`
@@ -94,9 +97,8 @@ impl Vary for () {
     }
     fn dv_dt(&self, _: &Self, _: f32) {}
     fn step(&self, _: &Self::Diff) {}
-
-    fn z_div(&self, _: f32) {}
 }
+impl ZDiv for () {}
 
 impl<T: Vary, U: Vary> Vary for (T, U) {
     type Iter = Iter<Self>;
@@ -114,12 +116,19 @@ impl<T: Vary, U: Vary> Vary for (T, U) {
     fn step(&self, (d0, d1): &Self::Diff) -> Self {
         (self.0.step(d0), self.1.step(d1))
     }
-
-    fn z_div(&self, z: f32) -> Self {
+}
+impl<T: ZDiv, U: ZDiv> ZDiv for (T, U) {
+    fn z_div(self, z: f32) -> Self {
         (self.0.z_div(z), self.1.z_div(z))
     }
 }
 
+impl ZDiv for f32 {
+    fn z_div(self, z: f32) -> Self {
+        self / z
+    }
+}
+
 impl<T: Vary> Iterator for Iter<T> {
     type Item = T;
 
diff --git a/core/src/math/vec.rs b/core/src/math/vec.rs
index c5432ebc..ffda8b8d 100644
--- a/core/src/math/vec.rs
+++ b/core/src/math/vec.rs
@@ -5,13 +5,15 @@
 use core::array;
 use core::fmt::{Debug, Formatter};
 use core::iter::Sum;
-use core::marker::PhantomData;
+use core::marker::PhantomData as Pd;
 use core::ops::{Add, Div, Index, IndexMut, Mul, Neg, Sub};
 use core::ops::{AddAssign, DivAssign, MulAssign, SubAssign};
 
 use crate::math::approx::ApproxEq;
 use crate::math::float::f32;
+use crate::math::point::Point;
 use crate::math::space::{Affine, Linear, Proj4, Real};
+use crate::math::vary::ZDiv;
 
 //
 // Types
@@ -30,11 +32,11 @@ use crate::math::space::{Affine, Linear, Proj4, Real};
 /// # Examples
 /// TODO
 #[repr(transparent)]
-pub struct Vector<Repr, Space = ()>(pub Repr, PhantomData<Space>);
+pub struct Vector<Repr, Space = ()>(pub Repr, Pd<Space>);
 
 /// A 2-vector with `f32` components.
 pub type Vec2<Basis = ()> = Vector<[f32; 2], Real<2, Basis>>;
-/// A 2-vector with `f32` components.
+/// A 3-vector with `f32` components.
 pub type Vec3<Basis = ()> = Vector<[f32; 3], Real<3, Basis>>;
 /// A `f32` 4-vector in the projective 3-space over ℝ, aka P<sub>3</sub>(ℝ).
 pub type ProjVec4 = Vector<[f32; 4], Proj4>;
@@ -54,12 +56,12 @@ pub type Vec2u<Basis = ()> = Vector<[u32; 2], Real<2, Basis>>;
 
 /// Returns a real 2-vector with components `x` and `y`.
 pub const fn vec2<Sc, B>(x: Sc, y: Sc) -> Vector<[Sc; 2], Real<2, B>> {
-    Vector([x, y], PhantomData)
+    Vector([x, y], Pd)
 }
 
 /// Returns a real 3-vector with components `x`, `y`, and `z`.
 pub const fn vec3<Sc, B>(x: Sc, y: Sc, z: Sc) -> Vector<[Sc; 3], Real<3, B>> {
-    Vector([x, y, z], PhantomData)
+    Vector([x, y, z], Pd)
 }
 
 /// Returns a vector with all components equal to a scalar.
@@ -74,7 +76,7 @@ pub const fn vec3<Sc, B>(x: Sc, y: Sc, z: Sc) -> Vector<[Sc; 3], Real<3, B>> {
 /// assert_eq!(v, vec3(1.23, 1.23, 1.23));
 #[inline]
 pub fn splat<Sp, Sc: Clone, const DIM: usize>(s: Sc) -> Vector<[Sc; DIM], Sp> {
-    s.into()
+    array::from_fn(|_| s.clone()).into()
 }
 
 //
@@ -85,7 +87,7 @@ impl<R, Sp> Vector<R, Sp> {
     /// Returns a new vector with representation `repr`.
     #[inline]
     pub const fn new(repr: R) -> Self {
-        Self(repr, PhantomData)
+        Self(repr, Pd)
     }
 
     /// Returns a vector with value equal to `self` but in space `S`.
@@ -99,6 +101,13 @@ impl<R, Sp> Vector<R, Sp> {
     pub fn to<S>(self) -> Vector<R, S> {
         Vector::new(self.0)
     }
+
+    /// Returns the affine point equivalent to `self`.
+    // TODO Cannot be const (yet?) due to E0493 :(
+    #[inline]
+    pub fn to_pt(self) -> Point<R, Sp> {
+        Point::new(self.0)
+    }
 }
 
 // TODO Many of these functions could be more generic
@@ -147,6 +156,9 @@ impl<Sp, const N: usize> Vector<[f32; N], Sp> {
     /// // Clamp to the unit cube
     /// let v = v.clamp(&splat(-1.0), &splat(1.0));
     /// assert_eq!(v, vec3(0.5, 1.0, -1.0));
+    // TODO f32 and f64 have inherent clamp methods because they're not Ord.
+    //      A generic clamp for Sc: Ord would conflict with this one. There is
+    //      currently no clean way to support both floats and impl Ord types.
     #[must_use]
     pub fn clamp(&self, min: &Self, max: &Self) -> Self {
         array::from_fn(|i| self[i].clamp(min[i], max[i])).into()
@@ -206,13 +218,15 @@ where
     {
         other.mul(self.scalar_project(other))
     }
+}
 
+impl<Sc: Copy, Sp, const N: usize> Vector<[Sc; N], Sp> {
     /// Returns a vector of the same dimension as `self` by applying `f`
     /// component-wise.
     #[inline]
     #[must_use]
     pub fn map<T>(self, mut f: impl FnMut(Sc) -> T) -> Vector<[T; N], Sp> {
-        array::from_fn(|i| f(self[i])).into()
+        array::from_fn(|i| f(self.0[i])).into()
     }
 }
 
@@ -325,10 +339,10 @@ where
 
 impl<Sc, Sp, const DIM: usize> Affine for Vector<[Sc; DIM], Sp>
 where
-    Sc: Affine,
-    Sc::Diff: Linear<Scalar = Sc::Diff> + Copy,
+    Sc: Affine<Diff: Linear<Scalar = Sc::Diff> + Copy>,
 {
     type Space = Sp;
+    // TODO Vectors always Linear once Point used for affine stuff
     type Diff = Vector<[Sc::Diff; DIM], Sp>;
 
     /// The dimension (number of components) of `Self`.
@@ -337,17 +351,16 @@ where
     #[inline]
     fn add(&self, other: &Self::Diff) -> Self {
         // TODO Profile performance of array::from_fn
-        array::from_fn(|i| self.0[i].add(&other.0[i])).into()
+        Self(array::from_fn(|i| self.0[i].add(&other.0[i])), Pd)
     }
     #[inline]
     fn sub(&self, other: &Self) -> Self::Diff {
-        array::from_fn(|i| self.0[i].sub(&other.0[i])).into()
+        Vector(array::from_fn(|i| self.0[i].sub(&other.0[i])), Pd)
     }
 }
 
 impl<Sc, Sp, const DIM: usize> Linear for Vector<[Sc; DIM], Sp>
 where
-    Self: Affine<Diff = Self>,
     Sc: Linear<Scalar = Sc> + Copy,
 {
     type Scalar = Sc;
@@ -359,11 +372,20 @@ where
     }
     #[inline]
     fn neg(&self) -> Self {
-        array::from_fn(|i| self.0[i].neg()).into()
+        self.map(|c| c.neg())
     }
     #[inline]
     fn mul(&self, scalar: Self::Scalar) -> Self {
-        array::from_fn(|i| self.0[i].mul(scalar)).into()
+        self.map(|c| c.mul(scalar))
+    }
+}
+
+impl<Sc, Sp, const N: usize> ZDiv for Vector<[Sc; N], Sp>
+where
+    Sc: ZDiv + Copy,
+{
+    fn z_div(self, z: f32) -> Self {
+        self.map(|c| c.z_div(z))
     }
 }
 
@@ -389,13 +411,13 @@ impl<R: Copy, S> Copy for Vector<R, S> {}
 
 impl<R: Clone, S> Clone for Vector<R, S> {
     fn clone(&self) -> Self {
-        Self(self.0.clone(), PhantomData)
+        Self(self.0.clone(), Pd)
     }
 }
 
 impl<R: Default, S> Default for Vector<R, S> {
     fn default() -> Self {
-        Self(R::default(), PhantomData)
+        Self(R::default(), Pd)
     }
 }
 
@@ -407,21 +429,6 @@ impl<R: PartialEq, S> PartialEq for Vector<R, S> {
     }
 }
 
-impl<const DIM: usize, B> Debug for Real<DIM, B>
-where
-    B: Debug + Default,
-{
-    fn fmt(&self, f: &mut Formatter<'_>) -> core::fmt::Result {
-        const DIMS: [char; 3] = ['²', '³', '⁴'];
-        let b = B::default();
-        if let Some(dim) = DIMS.get(DIM - 2) {
-            write!(f, "ℝ{dim}<{b:?}>")
-        } else {
-            write!(f, "ℝ^{DIM}<{b:?}>")
-        }
-    }
-}
-
 impl<R: Debug, Sp: Debug + Default> Debug for Vector<R, Sp> {
     fn fmt(&self, f: &mut Formatter<'_>) -> core::fmt::Result {
         write!(f, "Vec<{:?}>", Sp::default())?;
@@ -431,8 +438,8 @@ impl<R: Debug, Sp: Debug + Default> Debug for Vector<R, Sp> {
 
 impl<R, Sp> From<R> for Vector<R, Sp> {
     #[inline]
-    fn from(els: R) -> Self {
-        Self(els, PhantomData)
+    fn from(repr: R) -> Self {
+        Self(repr, Pd)
     }
 }
 
@@ -442,7 +449,7 @@ impl<Sp, Sc: Clone, const DIM: usize> From<Sc> for Vector<[Sc; DIM], Sp> {
     /// This operation is also called "splat" or "broadcast".
     #[inline]
     fn from(scalar: Sc) -> Self {
-        array::from_fn(|_| scalar.clone()).into()
+        splat(scalar)
     }
 }
 
@@ -562,7 +569,7 @@ where
 {
     #[inline]
     fn div_assign(&mut self, rhs: f32) {
-        debug_assert!(f32::abs(rhs) > 1e-7);
+        debug_assert!(f32::abs(rhs) > 1e-7, "divisor {rhs} < epsilon");
         *self = Linear::mul(&*self, rhs.recip());
     }
 }
diff --git a/core/src/render.rs b/core/src/render.rs
index ffc14e3e..475e1aee 100644
--- a/core/src/render.rs
+++ b/core/src/render.rs
@@ -12,12 +12,13 @@ use crate::geom::{Tri, Vertex};
 use crate::math::{
     mat::{Mat4x4, RealToProj, RealToReal},
     vary::Vary,
-    vec::{vec3, ProjVec4, Vec3},
+    vec::{vec3, ProjVec4},
+    Lerp,
 };
 
-use clip::{view_frustum, Clip, ClipVert};
+use clip::{view_frustum, ClipVert};
 use ctx::{Context, DepthSort, FaceCull};
-use raster::tri_fill;
+use raster::{tri_fill, ScreenPt};
 use shader::{FragmentShader, VertexShader};
 use stats::Stats;
 use target::Target;
@@ -82,7 +83,7 @@ impl<S, Vtx, Var, Uni> Shader<Vtx, Var, Uni> for S where
 }
 
 /// Renders the given triangles into `target`.
-pub fn render<Vtx: Clone, Var: Vary, Uni: Copy, Shd>(
+pub fn render<Vtx: Clone, Var: Lerp + Vary, Uni: Copy, Shd>(
     tris: impl AsRef<[Tri<usize>]>,
     verts: impl AsRef<[Vtx]>,
     shader: &Shd,
@@ -93,32 +94,34 @@ pub fn render<Vtx: Clone, Var: Vary, Uni: Copy, Shd>(
 ) where
     Shd: Shader<Vtx, Var, Uni>,
 {
+    let verts = verts.as_ref();
+    let tris = tris.as_ref();
     let mut stats = Stats::start();
 
     stats.calls = 1.0;
-    stats.prims.i += tris.as_ref().len();
-    stats.verts.i += verts.as_ref().len();
+    stats.prims.i += tris.len();
+    stats.verts.i += verts.len();
 
     // Vertex shader: transform vertices to clip space
     let verts: Vec<_> = verts
-        .as_ref()
         .iter()
         // TODO Pass vertex as ref to shader
         .cloned()
-        .map(|v| ClipVert::new(shader.shade_vertex(v, uniform)))
+        .map(|v| shader.shade_vertex(v, uniform))
+        .map(ClipVert::new)
         .collect();
 
     // Map triangle vertex indices to actual vertices
     let tris: Vec<_> = tris
-        .as_ref()
         .iter()
         .map(|Tri(vs)| Tri(vs.map(|i| verts[i].clone())))
         .collect();
 
     // Clip against the view frustum
     let mut clipped = vec![];
-    tris.clip(&view_frustum::PLANES, &mut clipped);
+    view_frustum::clip(&tris[..], &mut clipped);
 
+    // Optional depth sorting for use case such as transparency
     if let Some(d) = ctx.depth_sort {
         depth_sort(&mut clipped, d);
     }
@@ -133,16 +136,18 @@ pub fn render<Vtx: Clone, Var: Vary, Uni: Copy, Shd>(
             // of the original view-space depth. The interpolated reciprocal
             // is used in fragment processing for depth testing (larger values
             // are closer) and for perspective correction of the varyings.
-            let pos = vec3(x, y, 1.0).z_div(w);
+            let pos = vec3(x, y, 1.0) / w; // TODO
             Vertex {
                 // Viewport transform
-                pos: to_screen.apply(&pos),
+                pos: to_screen.apply(&pos).to_pt(),
                 // Perspective correction
                 attrib: v.attrib.z_div(w),
             }
         });
 
         // Back/frontface culling
+        //
+        // TODO This could also be done earlier, before or as part of clipping
         match ctx.face_cull {
             Some(FaceCull::Back) if is_backface(&vs) => continue,
             Some(FaceCull::Front) if !is_backface(&vs) => continue,
@@ -174,7 +179,7 @@ fn depth_sort<A>(tris: &mut [Tri<ClipVert<A>>], d: DepthSort) {
     });
 }
 
-fn is_backface<V>(vs: &[Vertex<Vec3<Screen>, V>]) -> bool {
+fn is_backface<V>(vs: &[Vertex<ScreenPt, V>]) -> bool {
     let v = vs[1].pos - vs[0].pos;
     let u = vs[2].pos - vs[0].pos;
     v[0] * u[1] - v[1] * u[0] > 0.0
diff --git a/core/src/render/batch.rs b/core/src/render/batch.rs
index 9de94c1e..c98efb76 100644
--- a/core/src/render/batch.rs
+++ b/core/src/render/batch.rs
@@ -3,8 +3,8 @@
 use alloc::vec::Vec;
 use core::borrow::Borrow;
 
-use crate::geom::{mesh, Mesh, Tri};
-use crate::math::{mat::Mat4x4, vary::Vary};
+use crate::geom::{Mesh, Tri, Vertex3};
+use crate::math::{mat::Mat4x4, vary::Vary, Lerp};
 
 use super::{ctx::Context, target::Target, NdcToScreen, Shader};
 
@@ -77,7 +77,7 @@ impl<Vtx, Uni, Shd, Tgt, Ctx> Batch<Vtx, Uni, Shd, Tgt, Ctx> {
     pub fn mesh<A: Clone>(
         self,
         mesh: &Mesh<A>,
-    ) -> Batch<mesh::Vertex<A>, Uni, Shd, Tgt, Ctx> {
+    ) -> Batch<Vertex3<A>, Uni, Shd, Tgt, Ctx> {
         let faces = mesh.faces.clone();
         let verts = mesh.verts.clone();
         update!(verts faces; self uniform shader viewport target ctx)
@@ -118,7 +118,7 @@ where
     Ctx: Borrow<Context>,
 {
     /// Renders this batch of geometry.
-    pub fn render<V: Vary>(&mut self)
+    pub fn render<V: Lerp + Vary>(&mut self)
     where
         Shd: Shader<Vtx, V, Uni>,
     {
diff --git a/core/src/render/cam.rs b/core/src/render/cam.rs
index 3f8382f6..07476031 100644
--- a/core/src/render/cam.rs
+++ b/core/src/render/cam.rs
@@ -9,9 +9,11 @@ use crate::math::{
     space::Linear,
     vary::Vary,
     vec::{vec2, Vec3},
+    Lerp,
 };
 use crate::util::{rect::Rect, Dims};
 
+use crate::math::point::Point3;
 #[cfg(feature = "fp")]
 use crate::math::{
     angle::Angle,
@@ -114,7 +116,7 @@ impl<M> Camera<M> {
     }
 
     /// Sets up orthographic projection.
-    pub fn orthographic(mut self, bounds: Range<Vec3>) -> Self {
+    pub fn orthographic(mut self, bounds: Range<Point3>) -> Self {
         self.project = orthographic(bounds.start, bounds.end);
         self
     }
@@ -127,7 +129,7 @@ impl<M: Mode> Camera<M> {
     }
 
     /// Renders the given geometry from the viewpoint of this camera.
-    pub fn render<B, Vtx: Clone, Var: Vary, Uni: Copy, Shd>(
+    pub fn render<B, Vtx: Clone, Var: Lerp + Vary, Uni: Copy, Shd>(
         &self,
         tris: impl AsRef<[Tri<usize>]>,
         verts: impl AsRef<[Vtx]>,
diff --git a/core/src/render/clip.rs b/core/src/render/clip.rs
index b204cf50..0cc1bfe9 100644
--- a/core/src/render/clip.rs
+++ b/core/src/render/clip.rs
@@ -13,15 +13,14 @@
 //! doing it for every scanline individually.
 //!
 
-use alloc::vec;
-use alloc::vec::Vec;
-
-use crate::geom::{Plane, Tri, Vertex};
-use crate::math::vary::Vary;
-use crate::math::vec::{ProjVec4, Vec3};
+use alloc::{vec, vec::Vec};
+use std::iter::zip;
 
 use view_frustum::{outcode, status};
 
+use crate::geom::{vertex, Tri, Vertex};
+use crate::math::{vec::ProjVec4, Lerp};
+
 /// Trait for types that can be [clipped][self] against planes.
 ///
 /// # Note to implementors
@@ -53,9 +52,6 @@ pub trait Clip {
 /// A vector in clip space.
 pub type ClipVec = ProjVec4;
 
-/// A plane in clip space.
-pub type ClipPlane = Plane<ClipVec>;
-
 /// A vertex in clip space.
 #[derive(Copy, Clone, Debug, PartialEq)]
 pub struct ClipVert<A> {
@@ -80,14 +76,17 @@ enum Status {
     Hidden,
 }
 
+impl Outcode {
+    const INSIDE: u8 = 0b00_11_11_11;
+}
+
+#[derive(Debug, Copy, Clone)]
+pub struct ClipPlane(ClipVec, u8);
+
 impl ClipPlane {
     /// Creates a clip plane given a normal and offset from origin.
-    ///
-    /// TODO Floating-point arithmetic is not permitted in const functions
-    ///      so the offset must be negated for now.
-    pub const fn new(normal: Vec3, neg_offset: f32) -> Self {
-        let [x, y, z] = normal.0;
-        Self(ClipVec::new([x, y, z, neg_offset]))
+    pub const fn new(x: f32, y: f32, z: f32, off: f32, bit: u8) -> Self {
+        Self(ClipVec::new([x, y, z, -off]), bit)
     }
 
     /// Returns the signed distance between `pt` and `self`.
@@ -112,6 +111,22 @@ impl ClipPlane {
         self.0.dot(pt)
     }
 
+    /// Computes the outcode bit for `pt`.
+    ///
+    /// The result is nonzero if `pt` is inside this plane, zero otherwise.
+    #[inline]
+    pub fn outcode(&self, pt: &ClipVec) -> u8 {
+        (self.signed_dist(pt) <= 0.0) as u8 * self.1
+    }
+
+    /// Checks the outcode of `v` against `self`.
+    ///
+    /// Returns `true` if this plane's outcode bit is set, `false` otherwise.
+    #[inline]
+    pub fn is_inside<A>(&self, v: &ClipVert<A>) -> bool {
+        self.1 & v.outcode.0 != 0
+    }
+
     /// Clips the convex polygon given by `verts_in` against `self` and
     /// returns the resulting vertices in `verts_out`.
     ///
@@ -129,21 +144,19 @@ impl ClipPlane {
     ///         `--__  \
     ///              `--c
     /// ```
-    pub fn clip_simple_polygon<A: Vary>(
+    pub fn clip_simple_polygon<A: Lerp + Clone>(
         &self,
         verts_in: &[ClipVert<A>],
         verts_out: &mut Vec<ClipVert<A>>,
     ) {
-        let mut verts = verts_in
-            .iter()
-            .chain(&verts_in[..1])
-            .map(|v| (v, self.signed_dist(&v.pos)));
+        let mut verts = verts_in.iter().chain(&verts_in[..1]);
 
-        let Some((mut v0, mut d0)) = &verts.next() else {
+        let Some(mut v0) = verts.next() else {
             return;
         };
-        for (v1, d1) in verts {
-            if d0 <= 0.0 {
+
+        for v1 in verts {
+            if self.is_inside(v0) {
                 // v0 is inside; emit it as-is. If v1 is also inside, we don't
                 // have to do anything; it is emitted on the next iteration.
                 verts_out.push((*v0).clone());
@@ -151,23 +164,29 @@ impl ClipPlane {
                 // v0 is outside, discard it. If v1 is also outside, we don't
                 // have to do anything; it is discarded on the next iteration.
             }
+            // TODO Doesn't use is_inside because it can't distinguish the case
+            //   where a vertex lies exactly on the plane. Though that's mostly
+            //   a theoretical edge case (heh).
+            let d0 = self.signed_dist(&v0.pos);
+            let d1 = self.signed_dist(&v1.pos);
             if d0 * d1 < 0.0 {
-                // Edge crosses the plane surface. Split the edge in two by
-                // interpolating and emitting a new vertex at intersection.
+                // The edge crosses the plane surface. Split the edge in two
+                // by interpolating and emitting a new vertex at intersection.
                 // The new vertex becomes one of the endpoints of a new "clip"
                 // edge coincident with the plane.
+                let d0 = self.signed_dist(&v0.pos);
+                let d1 = self.signed_dist(&v1.pos);
 
                 // `t` is the fractional distance from `v0` to the intersection
                 // point. If condition guarantees that `d1 - d0` is nonzero.
                 let t = -d0 / (d1 - d0);
 
-                verts_out.push(ClipVert {
-                    pos: v0.pos.lerp(&v1.pos, t),
-                    attrib: v0.attrib.lerp(&v1.attrib, t),
-                    outcode: Outcode(0b111111), // inside!
-                });
+                verts_out.push(ClipVert::new(vertex(
+                    v0.pos.lerp(&v1.pos, t),
+                    v0.attrib.lerp(&v1.attrib, t),
+                )));
             }
-            (v0, d0) = (v1, d1);
+            v0 = v1;
         }
     }
 }
@@ -186,30 +205,27 @@ impl ClipPlane {
 ///
 /// TODO Describe clip space
 pub mod view_frustum {
-    use crate::geom::Plane;
-    use crate::math::vec3;
-
     use super::*;
 
     /// The near, far, left, right, bottom, and top clipping planes,
     /// in that order.
     pub const PLANES: [ClipPlane; 6] = [
-        Plane::new(vec3(0.0, 0.0, -1.0), -1.0), // Near
-        Plane::new(vec3(0.0, 0.0, 1.0), -1.0),  // Far
-        Plane::new(vec3(-1.0, 0.0, 0.0), -1.0), // Left
-        Plane::new(vec3(1.0, 0.0, 0.0), -1.0),  // Right
-        Plane::new(vec3(0.0, -1.0, 0.0), -1.0), // Bottom
-        Plane::new(vec3(0.0, 1.0, 0.0), -1.0),  // Top
+        ClipPlane::new(0.0, 0.0, -1.0, 1.0, 1), // Near
+        ClipPlane::new(0.0, 0.0, 1.0, 1.0, 2),  // Far
+        ClipPlane::new(-1.0, 0.0, 0.0, 1.0, 4), // Left
+        ClipPlane::new(1.0, 0.0, 0.0, 1.0, 8),  // Right
+        ClipPlane::new(0.0, -1.0, 0.0, 1.0, 16), // Bottom
+        ClipPlane::new(0.0, 1.0, 0.0, 1.0, 32), // Top
     ];
 
+    /// Clips geometry against the standard view frustum.
+    pub fn clip<G: Clip + ?Sized>(geom: &G, out: &mut Vec<G::Item>) {
+        geom.clip(&PLANES, out);
+    }
+
     /// Returns the outcode of the given point.
     pub(super) fn outcode(pt: &ClipVec) -> Outcode {
-        // Top Btm Rgt Lft Far Near
-        //  1   2   4   8   16  32
-        let code = PLANES
-            .iter()
-            .fold(0, |code, p| code << 1 | (p.signed_dist(pt) <= 0.0) as u8);
-
+        let code = PLANES.iter().map(|p| p.outcode(pt)).sum();
         Outcode(code)
     }
 
@@ -218,11 +234,11 @@ pub mod view_frustum {
         let (all, any) = vs.iter().fold((!0, 0), |(all, any), v| {
             (all & v.outcode.0, any | v.outcode.0)
         });
-        if any != 0b111111 {
-            // If there's at least one plane that all vertices are outside of,
+        if any != Outcode::INSIDE {
+            // If there's at least one plane outside which all vertices are,
             // then the whole polygon is hidden
             Status::Hidden
-        } else if all == 0b111111 {
+        } else if all == Outcode::INSIDE {
             // If each vertex is inside all planes, the polygon is fully visible
             Status::Visible
         } else {
@@ -246,14 +262,14 @@ pub mod view_frustum {
 ///
 /// [^1]: Ivan Sutherland, Gary W. Hodgman: Reentrant Polygon Clipping.
 ///        Communications of the ACM, vol. 17, pp. 32–42, 1974
-pub fn clip_simple_polygon<'a, A: Vary>(
-    planes: &[Plane<ClipVec>],
+pub fn clip_simple_polygon<'a, A: Lerp + Clone>(
+    planes: &[ClipPlane],
     verts_in: &'a mut Vec<ClipVert<A>>,
     verts_out: &'a mut Vec<ClipVert<A>>,
 ) {
     debug_assert!(verts_out.is_empty());
 
-    for (i, p) in planes.iter().enumerate() {
+    for (p, i) in zip(planes, 0..) {
         p.clip_simple_polygon(verts_in, verts_out);
         if verts_out.is_empty() {
             // Nothing left to clip; the polygon was fully outside
@@ -267,16 +283,13 @@ pub fn clip_simple_polygon<'a, A: Vary>(
 }
 
 impl<V> ClipVert<V> {
-    pub fn new(v: Vertex<ClipVec, V>) -> Self {
-        ClipVert {
-            pos: v.pos,
-            attrib: v.attrib,
-            outcode: outcode(&v.pos),
-        }
+    pub fn new(Vertex { pos, attrib }: Vertex<ClipVec, V>) -> Self {
+        let outcode = outcode(&pos);
+        ClipVert { pos, attrib, outcode }
     }
 }
 
-impl<A: Vary> Clip for [Tri<ClipVert<A>>] {
+impl<A: Lerp + Clone> Clip for [Tri<ClipVert<A>>] {
     type Item = Tri<ClipVert<A>>;
 
     fn clip(&self, planes: &[ClipPlane], out: &mut Vec<Self::Item>) {
@@ -299,7 +312,7 @@ impl<A: Vary> Clip for [Tri<ClipVert<A>>] {
             verts_in.extend(vs.clone());
             clip_simple_polygon(planes, &mut verts_in, &mut verts_out);
 
-            if let Some((a, rest)) = verts_out.split_first() {
+            if let [a, rest @ ..] = &verts_out[..] {
                 // Clipping a triangle results in an n-gon, where n depends on
                 // how many planes the triangle intersects. Turn the resulting
                 // n-gon into a fan of triangles with common vertex `a`, for
@@ -328,6 +341,7 @@ impl<A: Vary> Clip for [Tri<ClipVert<A>>] {
 mod tests {
     use crate::geom::vertex;
     use crate::math::vary::Vary;
+    use std::dbg;
 
     use super::view_frustum::*;
     use super::*;
@@ -356,25 +370,26 @@ mod tests {
 
     #[test]
     fn outcode_inside() {
-        assert_eq!(outcode(&vec(0.0, 0.0, 0.0)).0, 0b111111);
-        assert_eq!(outcode(&vec(1.0, 0.0, 0.0)).0, 0b111111);
-        assert_eq!(outcode(&vec(0.0, -1.0, 0.0)).0, 0b111111);
-        assert_eq!(outcode(&vec(0.0, 1.0, 1.0)).0, 0b111111);
+        let inside = Outcode::INSIDE;
+        assert_eq!(outcode(&vec(0.0, 0.0, 0.0)).0, inside);
+        assert_eq!(outcode(&vec(1.0, 0.0, 0.0)).0, inside);
+        assert_eq!(outcode(&vec(0.0, -1.0, 0.0)).0, inside);
+        assert_eq!(outcode(&vec(0.0, 1.0, 1.0)).0, inside);
     }
 
     #[test]
     fn outcode_outside() {
         // Top Btm Rgt Lft Far Near
-        //  1   2   4   8   16  32
-
-        // Outside near == 32
-        assert_eq!(outcode(&vec(0.0, 0.0, -1.5)).0, 0b011111);
-        // Outside right == 4
-        assert_eq!(outcode(&vec(2.0, 0.0, 0.0)).0, 0b111011);
-        // Outside bottom == 2
-        assert_eq!(outcode(&vec(0.0, -1.01, 0.0)).0, 0b111101);
-        // Outside far left == 16|8
-        assert_eq!(outcode(&vec(-2.0, 0.0, 2.0)).0, 0b100111);
+        //  32  16  8   4   2   1
+
+        // Outside near == 1
+        assert_eq!(outcode(&vec(0.0, 0.0, -1.5)).0, 0b11_11_10);
+        // Outside right == 8
+        assert_eq!(outcode(&vec(2.0, 0.0, 0.0)).0, 0b11_01_11);
+        // Outside bottom == 16
+        assert_eq!(outcode(&vec(0.0, -1.01, 0.0)).0, 0b10_11_11);
+        // Outside far left == 2|4
+        assert_eq!(outcode(&vec(-2.0, 0.0, 2.0)).0, 0b11_10_01);
     }
 
     #[test]
@@ -509,7 +524,7 @@ mod tests {
     }
 
     #[test]
-    fn tri_clip_all_planes_fully_inside() {
+    fn tri_clip_against_frustum_fully_inside() {
         let tr = tri(
             vec(-1.0, -1.0, -1.0),
             vec(1.0, 1.0, 0.0),
@@ -520,7 +535,7 @@ mod tests {
         assert_eq!(res, &[tr]);
     }
     #[test]
-    fn tri_clip_all_planes_fully_outside() {
+    fn tri_clip_against_frustum_fully_outside() {
         //    z
         //    ^
         //    2-------0
@@ -539,7 +554,7 @@ mod tests {
         assert_eq!(res, &[]);
     }
     #[test]
-    fn tri_clip_all_planes_result_is_quad() {
+    fn tri_clip_against_frustum_result_is_quad() {
         //    z
         //    ^
         //    2
@@ -565,7 +580,7 @@ mod tests {
     }
 
     #[test]
-    fn tri_clip_all_planes_result_is_heptagon() {
+    fn tri_clip_against_frustum_result_is_heptagon() {
         //        z
         //        ^       2
         //        ·   /  /
@@ -588,9 +603,9 @@ mod tests {
     }
 
     #[test]
-    fn tri_clip_all_cases() {
+    fn tri_clip_against_frustum_all_cases() {
         // Methodically go through every possible combination of every
-        // vertex inside/outside of every plane, including degenerate cases.
+        // vertex inside/outside every plane, including degenerate cases.
 
         let xs = || (-2.0).vary(1.0, Some(5));
 
@@ -606,19 +621,32 @@ mod tests {
         });
 
         let mut in_tris = 0;
+        let mut in_degen = 0;
         let mut out_tris = [0; 8];
+        let mut out_degen = 0;
+        let mut out_total = 0;
         for tr in tris {
             let res = &mut vec![];
             [tr].clip(&PLANES, res);
             assert!(
                 res.iter().all(in_bounds),
-                "clip returned oob vertex:\n  from: {:#?}\n  to: {:#?}",
+                "clip returned oob or degenerate vertex:\n\
+                    input: {:#?}\n\
+                    output: {:#?}",
                 tr,
                 &res
             );
             in_tris += 1;
+            in_degen += is_degenerate(&tr) as u32;
             out_tris[res.len()] += 1;
+            out_total += res.len();
+            out_degen += res.iter().filter(|t| is_degenerate(t)).count()
         }
+        dbg!(in_tris);
+        dbg!(in_degen);
+        dbg!(out_degen);
+        dbg!(out_total);
+
         assert_eq!(in_tris, 5i32.pow(9));
         assert_eq!(
             out_tris,
@@ -626,9 +654,12 @@ mod tests {
         );
     }
 
-    fn in_bounds(tri: &Tri<ClipVert<f32>>) -> bool {
-        tri.0
-            .iter()
+    fn is_degenerate(Tri([a, b, c]): &Tri<ClipVert<f32>>) -> bool {
+        a.pos == b.pos || a.pos == c.pos || b.pos == c.pos
+    }
+
+    fn in_bounds(Tri(vs): &Tri<ClipVert<f32>>) -> bool {
+        vs.iter()
             .flat_map(|v| (v.pos / v.pos.w()).0)
             .all(|a| a.abs() <= 1.00001)
     }
diff --git a/core/src/render/raster.rs b/core/src/render/raster.rs
index b68a2874..4070ec45 100644
--- a/core/src/render/raster.rs
+++ b/core/src/render/raster.rs
@@ -15,14 +15,15 @@ use core::fmt::Debug;
 use core::ops::Range;
 
 use crate::geom::Vertex;
-use crate::math::{Vary, Vec3};
+use crate::math::point::Point3;
+use crate::math::Vary;
 
 use super::Screen;
 
 /// A fragment, or a single "pixel" in a rasterized primitive.
 #[derive(Clone, Debug)]
 pub struct Frag<V> {
-    pub pos: ScreenVec,
+    pub pos: ScreenPt,
     pub var: V,
 }
 
@@ -46,12 +47,12 @@ pub struct ScanlineIter<V: Vary> {
     n: u32,
 }
 
-/// Vector in screen space.
+/// Point in screen space.
 /// `x` and `y` are viewport pixel coordinates, `z` is depth.
-pub type ScreenVec = Vec3<Screen>;
+pub type ScreenPt = Point3<Screen>;
 
 /// Values to interpolate across a rasterized primitive.
-pub type Varyings<V> = (ScreenVec, V);
+pub type Varyings<V> = (ScreenPt, V);
 
 impl<V: Vary> Scanline<V> {
     pub fn fragments(&mut self) -> impl Iterator<Item = Frag<V>> + '_ {
@@ -106,7 +107,7 @@ impl<V: Vary> Iterator for ScanlineIter<V> {
 /// `scanline_fn` for each scanline. The scanlines are guaranteed to cover
 /// exactly those pixels whose center point lies inside the triangle. For more
 /// information on the scanline conversion, see [`scan`].
-pub fn tri_fill<V, F>(mut verts: [Vertex<ScreenVec, V>; 3], mut scanline_fn: F)
+pub fn tri_fill<V, F>(mut verts: [Vertex<ScreenPt, V>; 3], mut scanline_fn: F)
 where
     V: Vary,
     F: FnMut(Scanline<V>),
@@ -243,8 +244,8 @@ mod tests {
 
     use crate::assert_approx_eq;
     use crate::geom::vertex;
-    use crate::math::vary::Vary;
-    use crate::math::vec3;
+    use crate::math::point::pt3;
+    use crate::math::vary::{Vary, ZDiv};
     use crate::util::buf::Buf2;
 
     use super::{tri_fill, Frag, Scanline};
@@ -256,10 +257,10 @@ mod tests {
         let mut buf = Buf2::new((20, 10));
 
         let verts = [
-            vec3(8.0, 0.0, 0.0),
-            vec3(0.0, 6.0, 0.0),
-            vec3(14.0, 10.0, 0.0),
-            vec3(20.0, 3.0, 0.0),
+            pt3(8.0, 0.0, 0.0),
+            pt3(0.0, 6.0, 0.0),
+            pt3(14.0, 10.0, 0.0),
+            pt3(20.0, 3.0, 0.0),
         ]
         .map(|pos| vertex(pos, 0.0));
 
@@ -299,12 +300,8 @@ mod tests {
     #[test]
     fn gradient() {
         use core::fmt::Write;
-        let verts = [
-            vec3::<_, ()>(15.0, 2.0, 0.0),
-            vec3(2.0, 8.0, 1.0),
-            vec3(26.0, 14.0, 0.5),
-        ]
-        .map(|pos| vertex(vec3(pos.x(), pos.y(), 1.0), pos.z()));
+        let verts = [(15.0, 2.0, 0.0), (2.0, 8.0, 1.0), (26.0, 14.0, 0.5)]
+            .map(|(x, y, val)| vertex(pt3(x, y, 1.0), val));
 
         let expected = r"
               0
@@ -341,8 +338,8 @@ mod tests {
             y: 42,
             xs: 8..16,
             vs: Vary::vary_to(
-                (vec3(8.0, 42.0, 1.0 / w0), 3.0.z_div(w0)),
-                (vec3(16.0, 42.0, 1.0 / w1), 5.0.z_div(w1)),
+                (pt3(8.0, 42.0, 1.0 / w0), 3.0f32.z_div(w0)),
+                (pt3(16.0, 42.0, 1.0 / w1), 5.0f32.z_div(w1)),
                 8,
             ),
         };
@@ -359,7 +356,7 @@ mod tests {
         let mut x = 8.0;
 
         for ((Frag { pos, var }, z), v) in sl.fragments().zip(zs).zip(vars) {
-            assert_approx_eq!(pos, vec3(x, 42.0, z.recip()));
+            assert_approx_eq!(pos, pt3(x, 42.0, z.recip()));
             assert_approx_eq!(var, v);
 
             x += 1.0;
diff --git a/core/src/util/buf.rs b/core/src/util/buf.rs
index 3acc4070..48a5bcc9 100644
--- a/core/src/util/buf.rs
+++ b/core/src/util/buf.rs
@@ -3,14 +3,16 @@
 //! Useful for storing pixel data of any kind, among other things.
 
 use alloc::{vec, vec::Vec};
-use core::fmt::{self, Debug, Formatter};
-use core::iter;
-use core::ops::{Deref, DerefMut};
-
-use crate::util::Dims;
+use core::{
+    fmt::{self, Debug, Formatter},
+    iter,
+    ops::{Deref, DerefMut},
+};
 
 use inner::Inner;
 
+use crate::util::Dims;
+
 //
 // Traits
 //
@@ -43,16 +45,16 @@ pub trait AsMutSlice2<T> {
 ///
 /// # Examples
 /// ```
-/// # use retrofire_core::util::buf::*;
-/// # use retrofire_core::math::vec::*;
+/// # use retrofire_core::util::buf::Buf2;
+/// # use retrofire_core::math::point::pt2;
 /// // Elements initialized with `Default::default()`
 /// let mut buf = Buf2::new((4, 4));
-/// // Indexing with a 2D vector (x, y) yields element at row y, column x:
-/// buf[vec2(2, 1)] = 123;
+/// // Indexing with a 2D point (x, y) yields element at row y, column x:
+/// buf[pt2(2, 1)] = 123;
 /// // Indexing with an usize i yields row with index i as a slice:
-/// assert_eq!(buf[1usize], [0, 0, 123, 0]);
+/// assert_eq!(buf[1], [0, 0, 123, 0]);
 /// // Thus you can also do this, row first, column second:
-/// assert_eq!(buf[1usize][2], 123)
+/// assert_eq!(buf[1][2], 123)
 /// ```
 #[derive(Clone)]
 #[repr(transparent)]
@@ -203,8 +205,8 @@ impl<'a, T> Slice2<'a, T> {
     /// # use retrofire_core::util::buf::Slice2;
     /// let data = &[0, 1, 2, 3, 4, 5, 6];
     /// let slice = Slice2::new((2, 2), 3, data);
-    /// assert_eq!(&slice[0usize], &[0, 1]);
-    /// assert_eq!(&slice[1usize], &[3, 4]);
+    /// assert_eq!(&slice[0], &[0, 1]);
+    /// assert_eq!(&slice[1], &[3, 4]);
     /// ```
     /// Above, `slice` represents a 2×2 rectangle with stride 3, such that
     /// the first row maps to `data[0..2]` and the second to `data[3..5]`:
@@ -347,7 +349,10 @@ pub mod inner {
         ops::{Deref, DerefMut, Index, IndexMut, Range},
     };
 
-    use crate::{math::vec::Vec2u, util::rect::Rect, util::Dims};
+    use crate::{
+        math::point::Point2u,
+        util::{rect::Rect, Dims},
+    };
 
     use super::{AsSlice2, MutSlice2, Slice2};
 
@@ -506,7 +511,7 @@ pub mod inner {
 
         /// Returns a reference to the element at `pos`,
         /// or `None` if `pos` is out of bounds.
-        pub fn get(&self, pos: impl Into<Vec2u>) -> Option<&T> {
+        pub fn get(&self, pos: impl Into<Point2u>) -> Option<&T> {
             let [x, y] = pos.into().0;
             self.to_index_checked(x, y).map(|i| &self.data[i])
         }
@@ -607,7 +612,7 @@ pub mod inner {
 
         /// Returns a mutable reference to the element at `pos`,
         /// or `None` if `pos` is out of bounds.
-        pub fn get_mut(&mut self, pos: impl Into<Vec2u>) -> Option<&mut T> {
+        pub fn get_mut(&mut self, pos: impl Into<Point2u>) -> Option<&mut T> {
             let [x, y] = pos.into().0;
             self.to_index_checked(x, y)
                 .map(|i| &mut self.data[i])
@@ -662,7 +667,7 @@ pub mod inner {
     impl<T, D, Pos> Index<Pos> for Inner<T, D>
     where
         D: Deref<Target = [T]>,
-        Pos: Into<Vec2u>,
+        Pos: Into<Point2u>,
     {
         type Output = T;
 
@@ -680,7 +685,7 @@ pub mod inner {
     impl<T, D, Pos> IndexMut<Pos> for Inner<T, D>
     where
         D: DerefMut<Target = [T]>,
-        Pos: Into<Vec2u>,
+        Pos: Into<Point2u>,
     {
         /// Returns a mutable reference to the element at position `pos`.
         ///
@@ -734,8 +739,8 @@ mod tests {
         assert_eq!(buf[2usize], [2, 12, 22, 32]);
 
         assert_eq!(buf[[0, 0]], 0);
-        assert_eq!(buf[vec2(1, 0)], 10);
-        assert_eq!(buf[vec2(3, 4)], 34);
+        assert_eq!(buf[[1, 0]], 10);
+        assert_eq!(buf[[3, 4]], 34);
 
         assert_eq!(buf.get([2, 3]), Some(&23));
         assert_eq!(buf.get([4, 4]), None);
@@ -749,7 +754,7 @@ mod tests {
         buf[2usize][1] = 123;
         assert_eq!(buf[2usize], [2, 123, 22, 32]);
 
-        buf[vec2(2, 3)] = 234;
+        buf[[2, 3]] = 234;
         assert_eq!(buf[[2, 3]], 234);
 
         *buf.get_mut([3, 4]).unwrap() = 345;
@@ -952,12 +957,12 @@ mod tests {
         let buf = Buf2::new_with((5, 4), |x, y| x * 10 + y);
         let slice = buf.slice((2.., 1..3));
 
-        assert_eq!(slice[vec2(0, 0)], 21);
-        assert_eq!(slice[vec2(1, 0)], 31);
-        assert_eq!(slice[vec2(2, 1)], 42);
+        assert_eq!(slice[[0, 0]], 21);
+        assert_eq!(slice[[1, 0]], 31);
+        assert_eq!(slice[[2, 1]], 42);
 
-        assert_eq!(slice.get(vec2(2, 1)), Some(&42));
-        assert_eq!(slice.get(vec2(2, 2)), None);
+        assert_eq!(slice.get([2, 1]), Some(&42));
+        assert_eq!(slice.get([2, 2]), None);
     }
 
     #[test]
@@ -966,10 +971,10 @@ mod tests {
         let mut slice = buf.slice_mut((2.., 1..3));
 
         slice[[2, 1]] = 123;
-        assert_eq!(slice[vec2(2, 1)], 123);
+        assert_eq!(slice[[2, 1]], 123);
 
-        assert_eq!(slice.get_mut(vec2(2, 1)), Some(&mut 123));
-        assert_eq!(slice.get(vec2(2, 2)), None);
+        assert_eq!(slice.get_mut([2, 1]), Some(&mut 123));
+        assert_eq!(slice.get([2, 2]), None);
 
         buf[[2, 2]] = 321;
         let slice = buf.slice((1.., 2..));
diff --git a/demos/Cargo.toml b/demos/Cargo.toml
index e905f069..8902cdfb 100644
--- a/demos/Cargo.toml
+++ b/demos/Cargo.toml
@@ -40,7 +40,7 @@ required-features = ["minifb", "re-geom"]
 
 [[bin]]
 name = "crates"
-required-features = ["minifb", "re-geom"]
+required-features = ["sdl2", "re-geom"]
 
 [[bin]]
 name = "sprites"
diff --git a/demos/src/bin/bezier.rs b/demos/src/bin/bezier.rs
index d647c7d1..3cabc6bf 100644
--- a/demos/src/bin/bezier.rs
+++ b/demos/src/bin/bezier.rs
@@ -1,15 +1,17 @@
 use std::mem::swap;
 use std::ops::ControlFlow::Continue;
 
-use re::math::rand::{Distrib, Uniform, UnitDisk, Xorshift64};
-use re::math::spline::BezierSpline;
 use re::prelude::*;
-use re::util::Dims;
 
-use re_front::minifb::Window;
-use re_front::Frame;
+use re::math::{
+    point::{pt2, Point2, Point2u},
+    rand::{Distrib, Uniform, UnitDisk, Xorshift64},
+    spline::BezierSpline,
+};
 
-fn line([mut p0, mut p1]: [Vec2; 2]) -> impl Iterator<Item = Vec2u> {
+use re_front::{dims::SVGA_800_600, minifb::Window, Frame};
+
+fn line([mut p0, mut p1]: [Point2; 2]) -> impl Iterator<Item = Point2u> {
     if p0.y() > p1.y() {
         swap(&mut p0, &mut p1);
     }
@@ -22,12 +24,14 @@ fn line([mut p0, mut p1]: [Vec2; 2]) -> impl Iterator<Item = Vec2u> {
         (vec2(dx / dy, 1.0), dy)
     };
 
-    p0.vary(step, Some(n as u32))
-        .map(|p| vec2(p.x() as u32, p.y() as u32))
+    // TODO Fix once points are Vary
+    p0.to_vec()
+        .vary(step, Some(n as u32))
+        .map(|p| p.map(|c| c as u32).to_pt())
 }
 
 fn main() {
-    let dims @ (w, h) = (640, 480);
+    let dims @ (w, h) = SVGA_800_600;
 
     let mut win = Window::builder()
         .title("retrofire//bezier")
@@ -39,8 +43,11 @@ fn main() {
     let pos = Uniform(vec2(0.0, 0.0)..vec2(w as f32, h as f32));
     let vel = UnitDisk;
 
-    let (mut pts, mut deltas): (Vec<Vec2>, Vec<Vec2>) =
-        (pos, vel).samples(rng).take(4).unzip();
+    let (mut pts, mut deltas): (Vec<Point2>, Vec<Vec2>) = (pos, vel)
+        .samples(rng)
+        .take(4)
+        .map(|(p, v)| (p.to_pt(), v))
+        .unzip();
 
     win.run(|Frame { dt, buf, .. }| {
         let b = BezierSpline::new(&pts);
@@ -54,16 +61,16 @@ fn main() {
             }
         }
 
-        let max = vec2((w - 1) as f32, (h - 1) as f32);
+        let max = pt2((w - 1) as f32, (h - 1) as f32);
         let secs = dt.as_secs_f32();
         for (p, d) in pts.iter_mut().zip(deltas.iter_mut()) {
-            *p = (*p + *d * 200.0 * secs).clamp(&Vec2::zero(), &max);
-
-            if p[0] == 0.0 || p[0] == max.x() {
-                d[0] = -d[0];
+            *p = (*p + *d * 200.0 * secs).clamp(&pt2(0.0, 0.0), &max);
+            let [dx, dy] = &mut d.0;
+            if p.x() == 0.0 || p.x() == max.x() {
+                *dx = -*dx;
             }
-            if p[1] == 0.0 || p[1] == max.y() {
-                d[1] = -d[1];
+            if p.y() == 0.0 || p.y() == max.y() {
+                *dy = -*dy;
             }
         }
         Continue(())
diff --git a/demos/src/bin/crates.rs b/demos/src/bin/crates.rs
index a600cf6c..fce3e483 100644
--- a/demos/src/bin/crates.rs
+++ b/demos/src/bin/crates.rs
@@ -2,12 +2,14 @@ use core::ops::ControlFlow::*;
 
 use re::prelude::*;
 
-use re::math::color::gray;
+use re::geom::Vertex3;
+use re::math::{color::gray, point::pt3};
 use re::render::{
     batch::Batch,
     cam::{Camera, FirstPerson},
     ModelToProj,
 };
+
 use re_front::sdl2::Window;
 use re_geom::solids::*;
 
@@ -18,13 +20,13 @@ fn main() {
         .expect("should create window");
 
     let floor_shader = Shader::new(
-        |v: Vertex<_, _>, mvp: &Mat4x4<ModelToProj>| {
+        |v: Vertex3<_>, mvp: &Mat4x4<ModelToProj>| {
             vertex(mvp.apply(&v.pos), v.attrib)
         },
         |frag: Frag<Color3f>| frag.var.to_color4(),
     );
     let crate_shader = Shader::new(
-        |v: Vertex<_, _>, mvp: &Mat4x4<ModelToProj>| {
+        |v: Vertex3<_>, mvp: &Mat4x4<ModelToProj>| {
             vertex(mvp.apply(&v.pos), v.attrib)
         },
         |frag: Frag<Normal3>| {
@@ -119,7 +121,7 @@ fn floor() -> Mesh<Color3f> {
         for i in -size..=size {
             let even_odd = ((i & 1) ^ (j & 1)) == 1;
 
-            let pos = vec3(i as f32, -1.0, j as f32);
+            let pos = pt3(i as f32, -1.0, j as f32);
             let col = if even_odd { gray(0.2) } else { gray(0.9) };
             bld.push_vert(pos, col);
 
diff --git a/demos/src/bin/solids.rs b/demos/src/bin/solids.rs
index 4b347769..14ea9058 100644
--- a/demos/src/bin/solids.rs
+++ b/demos/src/bin/solids.rs
@@ -4,11 +4,13 @@ use minifb::{Key, KeyRepeat};
 
 use re::prelude::*;
 
+use re::geom::Vertex3;
 use re::math::{
-    color::gray, mat::RealToReal, spline::smootherstep, vec::ProjVec4,
+    color::gray, mat::RealToReal, point::pt2, spline::smootherstep,
+    vec::ProjVec4,
 };
-use re::render::batch::Batch;
-use re::render::{cam::Camera, Model, ModelToProj, ModelToWorld};
+use re::render::{batch::Batch, cam::Camera, ModelToProj, ModelToWorld};
+
 use re_front::{minifb::Window, Frame};
 use re_geom::{io::parse_obj, solids::*};
 
@@ -62,7 +64,7 @@ fn main() {
         .perspective(1.5, 0.1..1000.0)
         .viewport(vec2(10, 10)..vec2(w - 10, h - 10));
 
-    type VertexIn = Vertex<Vec3<Model>, Normal3>;
+    type VertexIn = Vertex3<Normal3>;
     type VertexOut = Vertex<ProjVec4, Color3f>;
     type Uniform<'a> = (&'a Mat4x4<ModelToProj>, &'a Mat4x4<RealToReal<3>>);
 
@@ -156,11 +158,11 @@ fn objects(res: u32) -> [Mesh<Normal3>; 13] {
 fn lathe(secs: u32) -> Mesh<Normal3> {
     Lathe::new(
         vec![
-            vertex(vec2(0.75, -0.5), vec2(1.0, 1.0)),
-            vertex(vec2(0.55, -0.25), vec2(1.0, 0.5)),
-            vertex(vec2(0.5, 0.0), vec2(1.0, 0.0)),
-            vertex(vec2(0.55, 0.25), vec2(1.0, -0.5)),
-            vertex(vec2(0.75, 0.5), vec2(1.0, 1.0)),
+            vertex(pt2(0.75, -0.5), vec2(1.0, 1.0)),
+            vertex(pt2(0.55, -0.25), vec2(1.0, 0.5)),
+            vertex(pt2(0.5, 0.0), vec2(1.0, 0.0)),
+            vertex(pt2(0.55, 0.25), vec2(1.0, -0.5)),
+            vertex(pt2(0.75, 0.5), vec2(1.0, 1.0)),
         ],
         secs,
     )
diff --git a/demos/src/bin/sprites.rs b/demos/src/bin/sprites.rs
index 1c8885c4..f446b209 100644
--- a/demos/src/bin/sprites.rs
+++ b/demos/src/bin/sprites.rs
@@ -1,5 +1,6 @@
 use core::array::from_fn;
 use core::ops::ControlFlow::Continue;
+use re::geom::Vertex3;
 
 use re::prelude::*;
 
@@ -19,10 +20,11 @@ fn main() {
     ];
     let count = 10000;
     let rng = Xorshift64::default();
-    let verts: Vec<Vertex<Vec3<_>, Vec2<_>>> = UnitBall
+    // TODO Points in unit ball and other distribs
+    let verts: Vec<Vertex3<Vec2<_>>> = UnitBall
         .samples(rng)
         .take(count)
-        .flat_map(|pos| verts.map(|v| vertex(pos.to(), v)))
+        .flat_map(|pos| verts.map(|v| vertex(pos.to().to_pt(), v)))
         .collect();
 
     let tris: Vec<_> = (0..count)
@@ -36,10 +38,10 @@ fn main() {
         .expect("should create window");
 
     let shader = Shader::new(
-        |v: Vertex<Vec3<_>, Vec2<_>>,
+        |v: Vertex3<Vec2<_>>,
          (mv, proj): (&Mat4x4<ModelToView>, &Mat4x4<ViewToProj>)| {
             let vertex_pos = 0.008 * vec3(v.attrib.x(), v.attrib.y(), 0.0);
-            let view_pos = mv.apply(&v.pos) + vertex_pos;
+            let view_pos = mv.apply_pt(&v.pos) + vertex_pos;
             vertex(proj.apply(&view_pos), v.attrib)
         },
         |frag: Frag<Vec2<_>>| {
diff --git a/demos/src/bin/square.rs b/demos/src/bin/square.rs
index 880c4c71..2f44cb32 100644
--- a/demos/src/bin/square.rs
+++ b/demos/src/bin/square.rs
@@ -2,15 +2,20 @@ use std::ops::ControlFlow::*;
 
 use re::prelude::*;
 
-use re::render::{ctx::Context, render, tex::SamplerClamp, Model, ModelToProj};
+use re::{
+    geom::Vertex3,
+    math::point::pt3,
+    render::{ctx::Context, render, tex::SamplerClamp, ModelToProj},
+};
+
 use re_front::minifb::Window;
 
 fn main() {
-    let verts: [Vertex<Vec3<Model>, TexCoord>; 4] = [
-        vertex(vec3(-1.0, -1.0, 0.0), uv(0.0, 0.0)),
-        vertex(vec3(-1.0, 1.0, 0.0), uv(0.0, 1.0)),
-        vertex(vec3(1.0, -1.0, 0.0), uv(1.0, 0.0)),
-        vertex(vec3(1.0, 1.0, 0.0), uv(1.0, 1.0)),
+    let verts: [Vertex3<TexCoord>; 4] = [
+        vertex(pt3(-1.0, -1.0, 0.0), uv(0.0, 0.0)),
+        vertex(pt3(-1.0, 1.0, 0.0), uv(0.0, 1.0)),
+        vertex(pt3(1.0, -1.0, 0.0), uv(1.0, 0.0)),
+        vertex(pt3(1.0, 1.0, 0.0), uv(1.0, 1.0)),
     ];
 
     let mut win = Window::builder()
@@ -32,7 +37,7 @@ fn main() {
     }));
 
     let shader = Shader::new(
-        |v: Vertex<_, _>, mvp: &Mat4x4<ModelToProj>| {
+        |v: Vertex3<_>, mvp: &Mat4x4<ModelToProj>| {
             vertex(mvp.apply(&v.pos), v.attrib)
         },
         |frag: Frag<_>| SamplerClamp.sample(&checker, frag.var),
diff --git a/demos/wasm/src/triangle.rs b/demos/wasm/src/triangle.rs
index 40c549bf..c8d62fe8 100644
--- a/demos/wasm/src/triangle.rs
+++ b/demos/wasm/src/triangle.rs
@@ -2,17 +2,17 @@ use core::ops::ControlFlow::*;
 
 use wasm_bindgen::prelude::*;
 
-use re::geom::{vertex, Tri, Vertex};
+use re::geom::{vertex, Tri, Vertex3};
 use re::math::{
     color::{rgba, Color4f},
     mat::{perspective, rotate_z, translate, viewport},
+    point::pt3,
     rads, vec2, vec3,
 };
 use re::render::{raster::Frag, render, shader::Shader, ModelToView};
 use re::util::Dims;
 
-use re_front::dims::SVGA_800_600;
-use re_front::wasm::Window;
+use re_front::{dims::SVGA_800_600, wasm::Window};
 
 // Entry point from JS
 #[wasm_bindgen(start)]
@@ -25,9 +25,9 @@ pub fn start() {
     win.ctx.color_clear = Some(rgba(0, 0, 0, 0x80));
 
     let vs = [
-        vertex(vec3(-2.0, 1.0, 0.0), rgba(1.0, 0.2, 0.1, 0.9)),
-        vertex(vec3(2.0, 2.0, 0.0), rgba(0.2, 0.9, 0.1, 0.8)),
-        vertex(vec3(0.0, -2.0, 0.0), rgba(0.3, 0.4, 1.0, 1.0)),
+        vertex(pt3(-2.0, 1.0, 0.0), rgba(1.0, 0.2, 0.1, 0.9)),
+        vertex(pt3(2.0, 2.0, 0.0), rgba(0.2, 0.9, 0.1, 0.8)),
+        vertex(pt3(0.0, -2.0, 0.0), rgba(0.3, 0.4, 1.0, 1.0)),
     ];
 
     let proj = perspective(1.0, 4.0 / 3.0, 0.1..1000.0);
@@ -42,7 +42,7 @@ pub fn start() {
         let mvp = mv.then(&proj);
 
         let sh = Shader::new(
-            |v: Vertex<_, Color4f>, _| vertex(mvp.apply(&v.pos), v.attrib),
+            |v: Vertex3<Color4f>, _| vertex(mvp.apply(&v.pos), v.attrib),
             |f: Frag<Color4f>| f.var.to_color4(),
         );
 
diff --git a/geom/src/io.rs b/geom/src/io.rs
index f16b52d0..0488dc26 100644
--- a/geom/src/io.rs
+++ b/geom/src/io.rs
@@ -37,11 +37,11 @@
 //! f 1 2 3 4 5
 //! ```
 
-use alloc::string::String;
-use alloc::vec;
-use core::fmt;
-use core::num::{ParseFloatError, ParseIntError};
-
+use alloc::{string::String, vec};
+use core::{
+    fmt,
+    num::{ParseFloatError, ParseIntError},
+};
 #[cfg(feature = "std")]
 use std::{
     error,
@@ -51,7 +51,10 @@ use std::{
 };
 
 use re::geom::{mesh::Builder, vertex, Mesh, Normal3, Tri};
-use re::math::vec::{vec3, Vec3, Vector};
+use re::math::{
+    point::Point3,
+    vec::{vec3, Vec3, Vector},
+};
 use re::render::{
     tex::{uv, TexCoord},
     Model,
@@ -143,7 +146,7 @@ pub fn parse_obj(src: impl IntoIterator<Item = u8>) -> Result<Builder<()>> {
         match item.as_bytes() {
             [b'#', ..] => continue, // Skip comment
 
-            b"v" => verts.push(parse_vector(tokens)?),
+            b"v" => verts.push(parse_point(tokens)?),
             b"vt" => tcrds.push(parse_texcoord(tokens)?),
             b"vn" => norms.push(parse_normal(tokens)?),
 
@@ -204,9 +207,8 @@ fn parse_texcoord<'a>(
     Ok(uv(u, v))
 }
 
-fn parse_vector<'a>(
-    i: &mut impl Iterator<Item = &'a str>,
-) -> Result<Vec3<Model>> {
+fn parse_vector<'a>(i: &mut impl Iterator<Item = &'a str>) -> Result<Vec3<Model>>
+{
     let x = next(i)?.parse()?;
     let y = next(i)?.parse()?;
     let z = next(i)?.parse()?;
@@ -217,6 +219,12 @@ fn parse_normal<'a>(i: &mut impl Iterator<Item = &'a str>) -> Result<Normal3> {
     parse_vector(i).map(Vector::to)
 }
 
+fn parse_point<'a>(
+    i: &mut impl Iterator<Item = &'a str>,
+) -> Result<Point3<Model>> {
+    parse_vector(i).map(Vector::to_pt)
+}
+
 fn parse_index(s: &str) -> Result<usize> {
     // OBJ has one-based indices
     Ok(s.parse::<usize>()? - 1)
@@ -291,8 +299,7 @@ impl From<ParseIntError> for Error {
 
 #[cfg(test)]
 mod tests {
-    use re::geom::Tri;
-    use re::math::vec3;
+    use re::{geom::Tri, math::point::pt3};
 
     use super::*;
 
@@ -312,7 +319,7 @@ v 1.0 2.0      0.0";
         let mesh = parse_obj(input).unwrap().build();
 
         assert_eq!(mesh.faces, vec![Tri([0, 1, 3]), Tri([3, 0, 2])]);
-        assert_eq!(mesh.verts[3].pos, vec3(1.0, 2.0, 0.0));
+        assert_eq!(mesh.verts[3].pos, pt3(1.0, 2.0, 0.0));
     }
 
     #[test]
diff --git a/geom/src/solids.rs b/geom/src/solids.rs
index eb3382a5..592639dc 100644
--- a/geom/src/solids.rs
+++ b/geom/src/solids.rs
@@ -1,16 +1,17 @@
 //! Mesh approximations of various geometric shapes.
 
-use core::array::from_fn;
-use core::ops::Range;
-
 use alloc::{vec, vec::Vec};
+use core::{array::from_fn, ops::Range};
 
-use re::geom::{mesh::Builder, vertex, Mesh, Normal2, Normal3, Vertex};
+use re::geom::{
+    mesh::Builder, vertex, Mesh, Normal2, Normal3, Vertex, Vertex2,
+};
 use re::math::{
     angle::{degs, polar, turns, Angle},
     mat::rotate_y,
+    point::{pt2, pt3, Point3},
     vary::Vary,
-    vec::{splat, vec2, vec3, Vec2, Vec3},
+    vec::{vec2, vec3, Vec3},
 };
 use re::render::tex::{uv, TexCoord};
 
@@ -47,9 +48,9 @@ pub struct Tetrahedron;
 #[derive(Copy, Clone, Debug)]
 pub struct Box {
     /// The left bottom near corner of the box.
-    pub left_bot_near: Vec3,
+    pub left_bot_near: Point3,
     /// The right top far corner of the box.
-    pub right_top_far: Vec3,
+    pub right_top_far: Point3,
 }
 
 /// Regular octahedron.
@@ -100,7 +101,7 @@ pub struct Icosahedron;
 #[derive(Clone, Debug, Default)]
 pub struct Lathe {
     /// The polyline defining the shape.
-    pub pts: Vec<Vertex<Vec2, Normal2>>,
+    pub pts: Vec<Vertex2<Normal2, ()>>,
     /// The number of facets used to approximate the surface of revolution.
     pub sectors: u32,
     /// Whether to add flat caps to both ends of the object. Has no effect
@@ -162,12 +163,12 @@ impl Tetrahedron {
         use re::math::float::f32;
         let sqrt = f32::sqrt;
         let coords = [
-            vec3(0.0, 1.0, 0.0),
-            vec3(sqrt(8.0 / 9.0), -1.0 / 3.0, 0.0),
-            vec3(-sqrt(2.0 / 9.0), -1.0 / 3.0, sqrt(2.0 / 3.0)),
-            vec3(-sqrt(2.0 / 9.0), -1.0 / 3.0, -sqrt(2.0 / 3.0)),
+            pt3(0.0, 1.0, 0.0),
+            pt3(sqrt(8.0 / 9.0), -1.0 / 3.0, 0.0),
+            pt3(-sqrt(2.0 / 9.0), -1.0 / 3.0, sqrt(2.0 / 3.0)),
+            pt3(-sqrt(2.0 / 9.0), -1.0 / 3.0, -sqrt(2.0 / 3.0)),
         ];
-        let norms = [-coords[3], -coords[1], -coords[2], -coords[0]];
+        let norms = [3, 1, 2, 0].map(|i| -coords[i].to_vec());
 
         let mut b = Mesh::builder();
 
@@ -238,9 +239,10 @@ impl Box {
 
     /// Returns a cube centered on the origin, with the given side length.
     pub fn cube(side_len: f32) -> Self {
+        let l = 0.5 * side_len;
         Self {
-            left_bot_near: splat(-0.5 * side_len),
-            right_top_far: splat(0.5 * side_len),
+            left_bot_near: pt3(-l, -l, -l),
+            right_top_far: pt3(l, l, l),
         }
     }
 
@@ -250,8 +252,8 @@ impl Box {
         b.push_faces(Self::FACES);
         for (pos_i, [norm_i, _uv_i]) in Self::VERTS {
             let pos = from_fn(|i| {
-                self.left_bot_near[i]
-                    .lerp(&self.right_top_far[i], Self::COORDS[pos_i][i])
+                self.left_bot_near.0[i]
+                    .lerp(&self.right_top_far.0[i], Self::COORDS[pos_i][i])
             });
             b.push_vert(pos.into(), Self::NORMS[norm_i]);
         }
@@ -260,13 +262,13 @@ impl Box {
 }
 
 impl Octahedron {
-    const COORDS: [Vec3; 6] = [
-        vec3(-1.0, 0.0, 0.0),
-        vec3(0.0, -1.0, 0.0),
-        vec3(0.0, 0.0, -1.0),
-        vec3(0.0, 1.0, 0.0),
-        vec3(0.0, 0.0, 1.0),
-        vec3(1.0, 0.0, 0.0),
+    const COORDS: [Point3; 6] = [
+        pt3(-1.0, 0.0, 0.0),
+        pt3(0.0, -1.0, 0.0),
+        pt3(0.0, 0.0, -1.0),
+        pt3(0.0, 1.0, 0.0),
+        pt3(0.0, 0.0, 1.0),
+        pt3(1.0, 0.0, 0.0),
     ];
     const NORMS: [Normal3; 8] = [
         vec3(-1.0, -1.0, -1.0),
@@ -377,7 +379,7 @@ impl Dodecahedron {
             b.push_face(i5, i5 + 2, i5 + 3);
             b.push_face(i5, i5 + 3, i5 + 4);
             for &j in face {
-                b.push_vert(Self::COORDS[j].normalize(), n);
+                b.push_vert(Self::COORDS[j].normalize().to_pt(), n);
             }
         }
         b.build()
@@ -422,7 +424,7 @@ impl Icosahedron {
             let n = Self::NORMALS[i].normalize();
             b.push_face(3 * i, 3 * i + 1, 3 * i + 2);
             for vi in *vs {
-                b.push_vert(Self::COORDS[vi].normalize(), n);
+                b.push_vert(Self::COORDS[vi].normalize().to_pt(), n);
             }
         }
         b.build()
@@ -430,7 +432,7 @@ impl Icosahedron {
 }
 
 impl Lathe {
-    pub fn new(pts: Vec<Vertex<Vec2, Normal2>>, sectors: u32) -> Self {
+    pub fn new(pts: Vec<Vertex2<Normal2, ()>>, sectors: u32) -> Self {
         assert!(sectors >= 3, "sectors must be at least 3, was {sectors}");
         Self {
             pts,
@@ -461,12 +463,12 @@ impl Lathe {
 
         // Create vertices
         for Vertex { pos, attrib: n } in &pts {
-            let mut pos = start.apply(&vec3(pos.x(), pos.y(), 0.0));
+            let mut pos = start.apply_pt(&pt3(pos.x(), pos.y(), 0.0));
             let mut norm = start.apply(&vec3(n.x(), n.y(), 0.0)).normalize();
 
             for _ in 0..=secs {
                 b.push_vert(pos, norm);
-                pos = rot.apply(&pos);
+                pos = rot.apply_pt(&pos);
                 norm = rot.apply(&norm);
             }
         }
@@ -521,7 +523,7 @@ impl Sphere {
         let pts = degs(-90.0)
             .vary_to(degs(90.0), segments)
             .map(|alt| polar(radius, alt).to_cart())
-            .map(|pos| vertex(pos, pos))
+            .map(|pos| vertex(pos.to_pt(), pos))
             .collect();
 
         Lathe::new(pts, sectors).build()
@@ -534,7 +536,7 @@ impl Torus {
         let pts = turns(0.0)
             .vary_to(turns(1.0), self.minor_sectors)
             .map(|alt| polar(self.minor_radius, alt).to_cart())
-            .map(|v| vertex(vec2(self.major_radius, 0.0) + v, v))
+            .map(|v| vertex(pt2(self.major_radius, 0.0) + v, v))
             .collect();
 
         Lathe::new(pts, self.major_sectors).build()
@@ -558,8 +560,8 @@ impl Cylinder {
 impl Cone {
     /// Builds the conical mesh.
     pub fn build(self) -> Mesh<Normal3> {
-        let base_pt = vec2(self.base_radius, -1.0);
-        let apex_pt = vec2(self.apex_radius, 1.0);
+        let base_pt = pt2(self.base_radius, -1.0);
+        let apex_pt = pt2(self.apex_radius, 1.0);
         let n = apex_pt - base_pt;
         let n = vec2(n.y(), -n.x());
         let pts = vec![vertex(base_pt, n), vertex(apex_pt, n)];
@@ -578,14 +580,14 @@ impl Capsule {
         let bottom_pts: Vec<_> = degs(-90.0)
             .vary_to(degs(0.0), cap_segments)
             .map(|alt| polar(radius, alt).to_cart())
-            .map(|v| vertex(vec2(0.0, -1.0) + v, v))
+            .map(|v| vertex(pt2(0.0, -1.0) + v, v))
             .collect();
 
         // Top hemisphere
         let top_pts = bottom_pts
             .iter()
             .map(|Vertex { pos, attrib: n }| {
-                vertex(vec2(pos.x(), -pos.y()), vec2(n.x(), -n.y()))
+                vertex(pt2(pos.x(), -pos.y()), vec2(n.x(), -n.y()))
             })
             .rev();