diff --git a/spatial/r2/box.go b/spatial/r2/box.go
new file mode 100644
index 00000000..6c988b3b
--- /dev/null
+++ b/spatial/r2/box.go
@@ -0,0 +1,111 @@
+// Copyright ©2022 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package r2
+
+import "math"
+
+// Box is a 2D bounding box. Well formed Boxes have
+// Min components smaller than Max components.
+type Box struct {
+	Min, Max Vec
+}
+
+// NewBox is shorthand for Box{Min:Vec{x0,y0}, Max:Vec{x1,y1}}.
+// The sides are swapped so that the resulting Box is well formed.
+func NewBox(x0, y0, x1, y1 float64) Box {
+	return Box{
+		Min: Vec{X: math.Min(x0, x1), Y: math.Min(y0, y1)},
+		Max: Vec{X: math.Max(x0, x1), Y: math.Max(y0, y1)},
+	}
+}
+
+// Size returns the size of the Box.
+func (a Box) Size() Vec {
+	return Sub(a.Max, a.Min)
+}
+
+// Center returns the center of the Box.
+func (a Box) Center() Vec {
+	return Scale(0.5, Add(a.Min, a.Max))
+}
+
+// IsEmpty returns true if a Box's volume is zero
+// or if a Min component is greater than its Max component.
+func (a Box) Empty() bool {
+	return a.Min.X >= a.Max.X || a.Min.Y >= a.Max.Y
+}
+
+// Vertices returns a slice of the 4 vertices
+// corresponding to each of the Box's corners.
+//
+// The order of the vertices is CCW in the XY plane starting at the box minimum.
+// If viewing box from +Z position the ordering is as follows:
+//  1. Bottom left.
+//  2. Bottom right.
+//  3. Top right.
+//  4. Top left.
+func (a Box) Vertices() []Vec {
+	return []Vec{
+		0: a.Min,
+		1: {a.Max.X, a.Min.Y},
+		2: a.Max,
+		3: {a.Min.X, a.Max.Y},
+	}
+}
+
+// Union returns a box enclosing both the receiver and argument Boxes.
+func (a Box) Union(b Box) Box {
+	if a.Empty() {
+		return b
+	}
+	if b.Empty() {
+		return a
+	}
+	return Box{
+		Min: minElem(a.Min, b.Min),
+		Max: maxElem(a.Max, b.Max),
+	}
+}
+
+// Add adds v to the bounding box components.
+// It is the equivalent of translating the Box by v.
+func (a Box) Add(v Vec) Box {
+	return Box{Add(a.Min, v), Add(a.Max, v)}
+}
+
+// Scale returns a new Box scaled by a size vector around its center.
+// The scaling is done element wise, which is to say the Box's X size is
+// scaled by v.X. Negative components of v are interpreted as zero.
+func (a Box) Scale(v Vec) Box {
+	v = maxElem(v, Vec{})
+	// TODO(soypat): Probably a better way to do this.
+	return centeredBox(a.Center(), mulElem(v, a.Size()))
+}
+
+// centeredBox creates a Box with a given center and size. Size's negative
+// components are interpreted as zero so that resulting box is well formed.
+func centeredBox(center, size Vec) Box {
+	size = maxElem(size, Vec{})
+	half := Scale(0.5, absElem(size))
+	return Box{Min: Sub(center, half), Max: Add(center, half)}
+}
+
+// Contains returns true if v is contained within the bounds of the Box.
+func (a Box) Contains(v Vec) bool {
+	if a.Empty() {
+		return v == a.Min && v == a.Max
+	}
+	return a.Min.X <= v.X && v.X <= a.Max.X &&
+		a.Min.Y <= v.Y && v.Y <= a.Max.Y
+}
+
+// Canon returns the canonical version of b. The returned Box has minimum
+// and maximum coordinates swapped if necessary so that it is well-formed.
+func (b Box) Canon() Box {
+	return Box{
+		Min: minElem(b.Min, b.Max),
+		Max: maxElem(b.Min, b.Max),
+	}
+}
diff --git a/spatial/r2/box_test.go b/spatial/r2/box_test.go
new file mode 100644
index 00000000..0e0ea869
--- /dev/null
+++ b/spatial/r2/box_test.go
@@ -0,0 +1,210 @@
+// Copyright ©2022 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package r2
+
+import (
+	"testing"
+
+	"golang.org/x/exp/rand"
+)
+
+func TestBoxContains(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 200; i++ {
+		b := randomBox(rnd)
+		for j := 0; j < 10; j++ {
+			contained := b.random(rnd)
+			if !b.Contains(contained) {
+				t.Error("bounding box should contain Vec")
+			}
+		}
+		uncontained := [4]Vec{
+			Add(b.Max, Vec{1, 0}),
+			Add(b.Max, Vec{0, 1}),
+			Sub(b.Min, Vec{1, 0}),
+			Sub(b.Min, Vec{0, 1}),
+		}
+		for _, unc := range uncontained {
+			if b.Contains(unc) {
+				t.Error("box should not contain vec")
+			}
+		}
+	}
+}
+
+func TestBoxUnion(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 200; i++ {
+		b1 := randomBox(rnd)
+		b2 := randomBox(rnd)
+		u := b1.Union(b2)
+		for j := 0; j < 10; j++ {
+			contained := b1.random(rnd)
+			if !u.Contains(contained) {
+				t.Error("union should contain b1's Vec")
+			}
+			contained = b2.random(rnd)
+			if !u.Contains(contained) {
+				t.Error("union should contain b2's Vec")
+			}
+		}
+		uncontained := [4]Vec{
+			Add(maxElem(b1.Max, b2.Max), Vec{1, 0}),
+			Add(maxElem(b1.Max, b2.Max), Vec{0, 1}),
+			Sub(minElem(b1.Min, b2.Min), Vec{1, 0}),
+			Sub(minElem(b1.Min, b2.Min), Vec{0, 1}),
+		}
+		for _, unc := range uncontained {
+			if !b1.Contains(unc) && !b2.Contains(unc) && u.Contains(unc) {
+				t.Error("union should not contain Vec")
+			}
+		}
+	}
+}
+
+func TestBoxCenter(t *testing.T) {
+	const tol = 1e-11
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 300; i++ {
+		b := randomBox(rnd)
+		center := b.Center()
+		size := b.Size()
+		newBox := centeredBox(center, size)
+		if b.Empty() {
+			t.Fatal("random box result must be well formed")
+		}
+		if !vecApproxEqual(b.Min, newBox.Min, tol) {
+			t.Errorf("min values of box not equal. got %g, expected %g", newBox.Min, b.Min)
+		}
+		if !vecApproxEqual(b.Max, newBox.Max, tol) {
+			t.Errorf("max values of box not equal. got %g, expected %g", newBox.Max, b.Max)
+		}
+	}
+}
+
+func TestBoxAdd(t *testing.T) {
+	const tol = 1e-14
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 12; i++ {
+		b := randomBox(rnd)
+		v := randomVec(rnd)
+		got := b.Add(v)
+		want := Box{Min: Add(b.Min, v), Max: Add(b.Max, v)}
+		if !vecApproxEqual(got.Min, want.Min, tol) {
+			t.Error("box min incorrect result")
+		}
+		if !vecApproxEqual(got.Max, want.Max, tol) {
+			t.Error("box max incorrect result")
+		}
+	}
+}
+
+func TestBoxScale(t *testing.T) {
+	const tol = 1e-11
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 300; i++ {
+		b := randomBox(rnd)
+		size := b.Size()
+
+		scaler := absElem(randomVec(rnd))
+		scaled := b.Scale(scaler)
+		gotScaler := divElem(scaled.Size(), size)
+		if !vecApproxEqual(scaler, gotScaler, tol) {
+			t.Errorf("got scaled %g, expected %g", gotScaler, scaler)
+		}
+		center := b.Center()
+		scaledCenter := scaled.Center()
+		if !vecApproxEqual(center, scaledCenter, tol) {
+			t.Error("scale modified center")
+		}
+	}
+}
+
+func TestBoxEmpty(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 300; i++ {
+		v := absElem(randomVec(rnd))
+		b := randomBox(rnd)
+		min := b.Min
+		max := b.Max
+		if !(Box{Min: min, Max: min}).Empty() {
+			t.Error("Box{min,min} should be empty")
+		}
+		if !(Box{Min: max, Max: max}).Empty() {
+			t.Error("Box{max,max} should be empty")
+		}
+		bmm := Box{Min: min, Max: Sub(min, v)}
+		if !bmm.Empty() {
+			t.Error("Box{min,min-v} should be empty")
+		} else if bmm.Canon().Empty() {
+			t.Error("Canonical box of Box{min,min-v} is not empty")
+		}
+		bMM := Box{Min: Add(max, v), Max: max}
+		if !bMM.Empty() {
+			t.Error("Box{max+v,max} should be empty")
+		} else if bmm.Canon().Empty() {
+			t.Error("Canonical box of Box{max+v,max} is not empty")
+		}
+	}
+}
+func TestBoxCanon(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 300; i++ {
+		b := randomBox(rnd)
+		badBox := Box{Min: b.Max, Max: b.Min}
+		canon := badBox.Canon()
+		if canon != b {
+			t.Error("swapped box canon should be equal to original box")
+		}
+	}
+}
+
+func TestBoxVertices(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 300; i++ {
+		b := randomBox(rnd)
+		gots := b.Vertices()
+		wants := goldenVertices(b)
+		if len(gots) != len(wants) {
+			t.Fatalf("bad length of vertices. expect 4, got %d", len(gots))
+		}
+		for j, want := range wants {
+			got := gots[j]
+			if !vecEqual(want, got) {
+				t.Errorf("%dth vertex not equal", j)
+			}
+		}
+	}
+}
+
+// randomBox returns a random valid bounding Box.
+func randomBox(rnd *rand.Rand) Box {
+	spatialScale := randomRange(0, 2000)
+	boxScale := randomRange(0.01, 1000)
+	return centeredBox(Scale(spatialScale, randomVec(rnd)), Scale(boxScale, absElem(randomVec(rnd))))
+}
+
+// Random returns a random point within the Box.
+// used to facilitate testing
+func (b Box) random(rnd *rand.Rand) Vec {
+	return Vec{
+		X: randomRange(b.Min.X, b.Max.X),
+		Y: randomRange(b.Min.Y, b.Max.Y),
+	}
+}
+
+// randomRange returns a random float64 [a,b)
+func randomRange(a, b float64) float64 {
+	return a + (b-a)*rand.Float64()
+}
+
+func goldenVertices(a Box) []Vec {
+	return []Vec{
+		0: a.Min,
+		1: {a.Max.X, a.Min.Y},
+		2: a.Max,
+		3: {a.Min.X, a.Max.Y},
+	}
+}
diff --git a/spatial/r2/vector.go b/spatial/r2/vector.go
index 34764948..2d06933d 100644
--- a/spatial/r2/vector.go
+++ b/spatial/r2/vector.go
@@ -76,11 +76,6 @@ func Cos(p, q Vec) float64 {
 	return Dot(p, q) / (Norm(p) * Norm(q))
 }
 
-// Box is a 2D bounding box.
-type Box struct {
-	Min, Max Vec
-}
-
 // Rotation describes a rotation in 2D.
 type Rotation struct {
 	sin, cos float64
@@ -112,3 +107,48 @@ func (r Rotation) Rotate(p Vec) Vec {
 func (r Rotation) isIdentity() bool {
 	return r.sin == 0 && r.cos == 1
 }
+
+// minElem returns a vector with the element-wise
+// minimum components of vectors a and b.
+func minElem(a, b Vec) Vec {
+	return Vec{
+		X: math.Min(a.X, b.X),
+		Y: math.Min(a.Y, b.Y),
+	}
+}
+
+// maxElem returns a vector with the element-wise
+// maximum components of vectors a and b.
+func maxElem(a, b Vec) Vec {
+	return Vec{
+		X: math.Max(a.X, b.X),
+		Y: math.Max(a.Y, b.Y),
+	}
+}
+
+// absElem returns the vector with components set to their absolute value.
+func absElem(a Vec) Vec {
+	return Vec{
+		X: math.Abs(a.X),
+		Y: math.Abs(a.Y),
+	}
+}
+
+// mulElem returns the Hadamard product between vectors a and b.
+//  v = {a.X*b.X, a.Y*b.Y, a.Z*b.Z}
+func mulElem(a, b Vec) Vec {
+	return Vec{
+		X: a.X * b.X,
+		Y: a.Y * b.Y,
+	}
+}
+
+// divElem returns the Hadamard product between vector a
+// and the inverse components of vector b.
+//  v = {a.X/b.X, a.Y/b.Y, a.Z/b.Z}
+func divElem(a, b Vec) Vec {
+	return Vec{
+		X: a.X / b.X,
+		Y: a.Y / b.Y,
+	}
+}
diff --git a/spatial/r2/vector_test.go b/spatial/r2/vector_test.go
index 18f9014b..05dbbf0f 100644
--- a/spatial/r2/vector_test.go
+++ b/spatial/r2/vector_test.go
@@ -8,6 +8,7 @@ import (
 	"math"
 	"testing"
 
+	"golang.org/x/exp/rand"
 	"gonum.org/v1/gonum/floats/scalar"
 )
 
@@ -170,6 +171,7 @@ func TestNorm2(t *testing.T) {
 }
 
 func TestUnit(t *testing.T) {
+	const tol = 1e-14
 	for _, test := range []struct {
 		v, want Vec
 	}{
@@ -185,7 +187,7 @@ func TestUnit(t *testing.T) {
 		{Vec{1e4, math.MaxFloat32 - 1}, Vec{0, 1}},
 	} {
 		got := Unit(test.v)
-		if !vecApproxEqual(got, test.want) {
+		if !vecApproxEqual(got, test.want, tol) {
 			t.Errorf(
 				"Unit(%v) = %v, want %v",
 				test.v, got, test.want,
@@ -201,6 +203,7 @@ func TestUnit(t *testing.T) {
 }
 
 func TestCos(t *testing.T) {
+	const tol = 1e-14
 	for _, test := range []struct {
 		v1, v2 Vec
 		want   float64
@@ -211,7 +214,6 @@ func TestCos(t *testing.T) {
 		{Vec{1, 0}, Vec{0, 1}, 0},
 		{Vec{1, 0}, Vec{-1, 0}, -1},
 	} {
-		tol := 1e-14
 		got := Cos(test.v1, test.v2)
 		if !scalar.EqualWithinAbs(got, test.want, tol) {
 			t.Errorf("cos(%v, %v)= %v, want %v",
@@ -222,6 +224,7 @@ func TestCos(t *testing.T) {
 }
 
 func TestRotate(t *testing.T) {
+	const tol = 1e-14
 	for _, test := range []struct {
 		v, q  Vec
 		alpha float64
@@ -237,7 +240,7 @@ func TestRotate(t *testing.T) {
 		{Vec{2, 2}, Vec{2, 0}, math.Pi, Vec{2, -2}},
 	} {
 		got := Rotate(test.v, test.alpha, test.q)
-		if !vecApproxEqual(got, test.want) {
+		if !vecApproxEqual(got, test.want, tol) {
 			t.Errorf(
 				"rotate(%v, %v, %v)= %v, want=%v",
 				test.v, test.alpha, test.q, got, test.want,
@@ -254,12 +257,26 @@ func vecIsNaNAny(v Vec) bool {
 	return math.IsNaN(v.X) || math.IsNaN(v.Y)
 }
 
-func vecApproxEqual(a, b Vec) bool {
-	const tol = 1e-14
+func vecApproxEqual(a, b Vec, tol float64) bool {
+	if tol == 0 {
+		return vecEqual(a, b)
+	}
 	if vecIsNaNAny(a) || vecIsNaNAny(b) {
 		return vecIsNaN(a) && vecIsNaN(b)
 	}
-
 	return scalar.EqualWithinAbs(a.X, b.X, tol) &&
 		scalar.EqualWithinAbs(a.Y, b.Y, tol)
 }
+
+func randomVec(rnd *rand.Rand) (v Vec) {
+	v.X = (rnd.Float64() - 0.5) * 20
+	v.Y = (rnd.Float64() - 0.5) * 20
+	return v
+}
+
+func vecEqual(a, b Vec) bool {
+	if vecIsNaNAny(a) || vecIsNaNAny(b) {
+		return vecIsNaN(a) && vecIsNaN(b)
+	}
+	return a == b
+}