diff --git a/spatial/r3/icosahedron_example_test.go b/spatial/r3/icosahedron_example_test.go
new file mode 100644
index 000000000..916e2cece
--- /dev/null
+++ b/spatial/r3/icosahedron_example_test.go
@@ -0,0 +1,104 @@
+// 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 r3_test
+
+import (
+	"fmt"
+
+	"gonum.org/v1/gonum/spatial/r3"
+)
+
+func ExampleTriangle_icosphere() {
+	// This example generates a 3D icosphere from
+	// a starting icosahedron by subdividing surfaces.
+	// See https://schneide.blog/2016/07/15/generating-an-icosphere-in-c/.
+	const subdivisions = 5
+	// vertices is a slice of r3.Vec
+	// triangles is a slice of [3]int indices
+	// referencing the vertices.
+	vertices, triangles := icosahedron()
+	for i := 0; i < subdivisions; i++ {
+		vertices, triangles = subdivide(vertices, triangles)
+	}
+	var faces []r3.Triangle
+	for _, t := range triangles {
+		var face r3.Triangle
+		for i := 0; i < 3; i++ {
+			face[i] = vertices[t[i]]
+		}
+		faces = append(faces, face)
+	}
+	fmt.Println(faces)
+	// The 3D rendering of the icosphere is left as an exercise to the reader.
+}
+
+// edgeIdx represents an edge of the icosahedron
+type edgeIdx [2]int
+
+func subdivide(vertices []r3.Vec, triangles [][3]int) ([]r3.Vec, [][3]int) {
+	// We generate a lookup table of all newly generated vertices so as to not
+	// duplicate new vertices. edgeIdx has lower index first.
+	lookup := make(map[edgeIdx]int)
+	var result [][3]int
+	for _, triangle := range triangles {
+		var mid [3]int
+		for edge := 0; edge < 3; edge++ {
+			lookup, mid[edge], vertices = subdivideEdge(lookup, vertices, triangle[edge], triangle[(edge+1)%3])
+		}
+		newTriangles := [][3]int{
+			{triangle[0], mid[0], mid[2]},
+			{triangle[1], mid[1], mid[0]},
+			{triangle[2], mid[2], mid[1]},
+			{mid[0], mid[1], mid[2]},
+		}
+		result = append(result, newTriangles...)
+	}
+	return vertices, result
+}
+
+// subdivideEdge takes the vertices list and indices first and second which
+// refer to the edge that will be subdivided.
+// lookup is a table of all newly generated vertices from
+// previous calls to subdivideEdge so as to not duplicate vertices.
+func subdivideEdge(lookup map[edgeIdx]int, vertices []r3.Vec, first, second int) (map[edgeIdx]int, int, []r3.Vec) {
+	key := edgeIdx{first, second}
+	if first > second {
+		// Swap to ensure edgeIdx always has lower index first.
+		key[0], key[1] = key[1], key[0]
+	}
+	vertIdx, vertExists := lookup[key]
+	if !vertExists {
+		// If edge not already subdivided add
+		// new dividing vertex to lookup table.
+		edge0 := vertices[first]
+		edge1 := vertices[second]
+		point := r3.Unit(r3.Add(edge0, edge1)) // vertex at a normalized position.
+		vertices = append(vertices, point)
+		vertIdx = len(vertices) - 1
+		lookup[key] = vertIdx
+	}
+	return lookup, vertIdx, vertices
+}
+
+// icosahedron returns an icosahedron mesh.
+func icosahedron() (vertices []r3.Vec, triangles [][3]int) {
+	const (
+		radiusSqrt = 1.0 // Example designed for unit sphere generation.
+		X          = radiusSqrt * .525731112119133606
+		Z          = radiusSqrt * .850650808352039932
+		N          = 0.0
+	)
+	return []r3.Vec{
+			{-X, N, Z}, {X, N, Z}, {-X, N, -Z}, {X, N, -Z},
+			{N, Z, X}, {N, Z, -X}, {N, -Z, X}, {N, -Z, -X},
+			{Z, X, N}, {-Z, X, N}, {Z, -X, N}, {-Z, -X, N},
+		}, [][3]int{
+			{0, 1, 4}, {0, 4, 9}, {9, 4, 5}, {4, 8, 5},
+			{4, 1, 8}, {8, 1, 10}, {8, 10, 3}, {5, 8, 3},
+			{5, 3, 2}, {2, 3, 7}, {7, 3, 10}, {7, 10, 6},
+			{7, 6, 11}, {11, 6, 0}, {0, 6, 1}, {6, 10, 1},
+			{9, 11, 0}, {9, 2, 11}, {9, 5, 2}, {7, 11, 2},
+		}
+}
diff --git a/spatial/r3/triangle.go b/spatial/r3/triangle.go
new file mode 100644
index 000000000..00ba26141
--- /dev/null
+++ b/spatial/r3/triangle.go
@@ -0,0 +1,117 @@
+// 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 r3
+
+import "math"
+
+// Triangle represents a triangle in 3D space and
+// is composed by 3 vectors corresponding to the position
+// of each of the vertices. Ordering of these vertices
+// decides the "normal" direction.
+// Inverting ordering of two vertices inverts the resulting direction.
+type Triangle [3]Vec
+
+// Centroid returns the intersection of the three medians of the triangle
+// as a point in space.
+func (t Triangle) Centroid() Vec {
+	return Scale(1.0/3.0, Add(Add(t[0], t[1]), t[2]))
+}
+
+// Normal returns the vector with direction
+// perpendicular to the Triangle's face and magnitude
+// twice that of the Triangle's area. The ordering
+// of the triangle vertices decides the normal's resulting
+// direction. The returned vector is not normalized.
+func (t Triangle) Normal() Vec {
+	s1, s2, _ := t.sides()
+	return Cross(s1, s2)
+}
+
+// IsDegenerate returns true if all of triangle's vertices are
+// within tol distance of its longest side.
+func (t Triangle) IsDegenerate(tol float64) bool {
+	longIdx := t.longIdx()
+	// calculate vertex distance from longest side
+	ln := line{t[longIdx], t[(longIdx+1)%3]}
+	dist := ln.distance(t[(longIdx+2)%3])
+	return dist <= tol
+}
+
+// longIdx returns index of the longest side. The sides
+// of the triangles are are as follows:
+//  - Side 0 formed by vertices 0 and 1
+//  - Side 1 formed by vertices 1 and 2
+//  - Side 2 formed by vertices 0 and 2
+func (t Triangle) longIdx() int {
+	sides := [3]Vec{Sub(t[1], t[0]), Sub(t[2], t[1]), Sub(t[0], t[2])}
+	len2 := [3]float64{Norm2(sides[0]), Norm2(sides[1]), Norm2(sides[2])}
+	longLen := len2[0]
+	longIdx := 0
+	if len2[1] > longLen {
+		longLen = len2[1]
+		longIdx = 1
+	}
+	if len2[2] > longLen {
+		longIdx = 2
+	}
+	return longIdx
+}
+
+// Area returns the surface area of the triangle.
+func (t Triangle) Area() float64 {
+	// Heron's Formula, see https://en.wikipedia.org/wiki/Heron%27s_formula.
+	// Also see William M. Kahan (24 March 2000). "Miscalculating Area and Angles of a Needle-like Triangle"
+	// for more discussion. https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf.
+	a, b, c := t.orderedLengths()
+	A := (c + (b + a)) * (a - (c - b))
+	A *= (a + (c - b)) * (c + (b - a))
+	return math.Sqrt(A) / 4
+}
+
+// orderedLengths returns the lengths of the sides of the triangle such that
+// a ≤ b ≤ c.
+func (t Triangle) orderedLengths() (a, b, c float64) {
+	s1, s2, s3 := t.sides()
+	l1 := Norm(s1)
+	l2 := Norm(s2)
+	l3 := Norm(s3)
+	// sort-3
+	if l2 < l1 {
+		l1, l2 = l2, l1
+	}
+	if l3 < l2 {
+		l2, l3 = l3, l2
+		if l2 < l1 {
+			l1, l2 = l2, l1
+		}
+	}
+	return l1, l2, l3
+}
+
+// sides returns vectors for each of the sides of t.
+func (t Triangle) sides() (Vec, Vec, Vec) {
+	return Sub(t[1], t[0]), Sub(t[2], t[1]), Sub(t[0], t[2])
+}
+
+// line is an infinite 3D line
+// defined by two points on the line.
+type line [2]Vec
+
+// vecOnLine takes a value between 0 and 1 to linearly
+// interpolate a point on the line.
+//  vecOnLine(0) returns l[0]
+//  vecOnLine(1) returns l[1]
+func (l line) vecOnLine(t float64) Vec {
+	lineDir := Sub(l[1], l[0])
+	return Add(l[0], Scale(t, lineDir))
+}
+
+// distance returns the minimum euclidean distance of point p
+// to the line.
+func (l line) distance(p Vec) float64 {
+	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+	num := Norm(Cross(Sub(p, l[0]), Sub(p, l[1])))
+	return num / Norm(Sub(l[1], l[0]))
+}
diff --git a/spatial/r3/triangle_test.go b/spatial/r3/triangle_test.go
new file mode 100644
index 000000000..689079cf4
--- /dev/null
+++ b/spatial/r3/triangle_test.go
@@ -0,0 +1,203 @@
+// 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 r3
+
+import (
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+)
+
+func TestTriangleDegenerate(t *testing.T) {
+	const (
+		// tol is how much closer the problematic
+		// vertex is placed to avoid floating point error
+		// for degeneracy calculation.
+		tol = 1e-12
+		// This is the argument to Degenerate and represents
+		// the minimum permissible distance between the triangle
+		// longest edge and the opposite vertex.
+		spatialTol = 1e-2
+	)
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 200; i++ {
+		// Generate a random line for the longest triangle side.
+		ln := line{randomVec(rnd), randomVec(rnd)}
+		lineDir := Sub(ln[1], ln[0])
+
+		perpendicular := Unit(Cross(lineDir, randomVec(rnd)))
+		// generate 3 permutations of needle triangles for
+		// each vertex. A needle triangle has two vertices
+		// very close to eachother an its third vertex far away.
+		var needle Triangle
+		for j := 0; j < 3; j++ {
+			needle[j] = ln[0]
+			needle[(j+1)%3] = ln[1]
+			needle[(j+2)%3] = Add(ln[1], Scale((1-tol)*spatialTol, perpendicular))
+			if !needle.IsDegenerate(spatialTol) {
+				t.Error("needle triangle not degenerate")
+			}
+		}
+
+		midpoint := ln.vecOnLine(0.5)
+		// cap triangles are characterized by having two sides
+		// of similar lengths and whose sum is approximately equal
+		// to the remaining longest side.
+		var cap Triangle
+		for j := 0; j < 3; j++ {
+			cap[j] = ln[0]
+			cap[(j+1)%3] = ln[1]
+			cap[(j+2)%3] = Add(midpoint, Scale((1-tol)*spatialTol, perpendicular))
+			if !cap.IsDegenerate(spatialTol) {
+				t.Error("cap triangle not degenerate")
+			}
+		}
+
+		var degenerate Triangle
+		for j := 0; j < 3; j++ {
+			degenerate[j] = ln[0]
+			degenerate[(j+1)%3] = ln[1]
+			// vertex perpendicular to some random point on longest side.
+			degenerate[(j+2)%3] = Add(ln.vecOnLine(rnd.Float64()), Scale((1-tol)*spatialTol, perpendicular))
+			if !degenerate.IsDegenerate(spatialTol) {
+				t.Error("random degenerate triangle not degenerate")
+			}
+			// vertex about longest side 0 vertex
+			degenerate[(j+2)%3] = Add(ln[0], Scale((1-tol)*spatialTol, Unit(randomVec(rnd))))
+			if !degenerate.IsDegenerate(spatialTol) {
+				t.Error("needle-like degenerate triangle not degenerate")
+			}
+			// vertex about longest side 1 vertex
+			degenerate[(j+2)%3] = Add(ln[1], Scale((1-tol)*spatialTol, Unit(randomVec(rnd))))
+			if !degenerate.IsDegenerate(spatialTol) {
+				t.Error("needle-like degenerate triangle not degenerate")
+			}
+		}
+	}
+}
+
+func TestTriangleCentroid(t *testing.T) {
+	const tol = 1e-12
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 100; i++ {
+		tri := randomTriangle(rnd)
+		got := tri.Centroid()
+		want := Vec{
+			X: (tri[0].X + tri[1].X + tri[2].X) / 3,
+			Y: (tri[0].Y + tri[1].Y + tri[2].Y) / 3,
+			Z: (tri[0].Z + tri[1].Z + tri[2].Z) / 3,
+		}
+		if !vecApproxEqual(got, want, tol) {
+			t.Fatalf("got %.6g, want %.6g", got, want)
+		}
+	}
+}
+
+func TestTriangleNormal(t *testing.T) {
+	const tol = 1e-12
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 100; i++ {
+		tri := randomTriangle(rnd)
+		got := tri.Normal()
+		expect := goldenNormal(tri)
+		if !vecApproxEqual(got, expect, tol) {
+			t.Fatalf("got %.6g, want %.6g", got, expect)
+		}
+	}
+}
+
+func TestTriangleArea(t *testing.T) {
+	const tol = 1e-16
+	for _, test := range []struct {
+		T      Triangle
+		Expect float64
+	}{
+		{
+			T: Triangle{
+				{0, 0, 0},
+				{1, 0, 0},
+				{0, 1, 0},
+			},
+			Expect: 0.5,
+		},
+		{
+			T: Triangle{
+				{1, 0, 0},
+				{0, 1, 0},
+				{0, 0, 0},
+			},
+			Expect: 0.5,
+		},
+		{
+			T: Triangle{
+				{20, 0, 0},
+				{0, 0, 20},
+				{0, 0, 0},
+			},
+			Expect: 20 * 20 / 2,
+		},
+	} {
+		got := test.T.Area()
+		if math.Abs(got-test.Expect) > tol {
+			t.Errorf("got area %g, expected %g", got, test.Expect)
+		}
+		if test.T.IsDegenerate(tol) {
+			t.Error("well-formed triangle is degenerate")
+		}
+	}
+	const tol2 = 1e-12
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 100; i++ {
+		tri := randomTriangle(rnd)
+		got := tri.Area()
+		want := Norm(Cross(Sub(tri[1], tri[0]), Sub(tri[2], tri[0]))) / 2
+		if math.Abs(got-want) > tol2 {
+			t.Errorf("got area %g not match half norm of cross product %g", got, want)
+		}
+	}
+}
+
+func TestTriangleOrderedLengths(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i := 0; i < 200; i++ {
+		tri := randomTriangle(rnd)
+		s1, s2, s3 := tri.sides()
+		l1 := Norm(s1)
+		l2 := Norm(s2)
+		l3 := Norm(s3)
+		a, b, c := tri.orderedLengths()
+		if a != l1 && a != l2 && a != l3 {
+			t.Error("shortest ordered length not a side of the triangle")
+		}
+		if b != l1 && b != l2 && b != l3 {
+			t.Error("middle ordered length not a side of the triangle")
+		}
+		if c != l1 && c != l2 && c != l3 {
+			t.Error("longest ordered length not a side of the triangle")
+		}
+		if a > b || a > c {
+			t.Error("ordered short side not shortest side")
+		}
+		if c < b {
+			t.Error("ordered long side not longest side")
+		}
+	}
+}
+
+// taken from soypat/sdf library where it has been thoroughly tested empirically.
+func goldenNormal(t Triangle) Vec {
+	e1 := Sub(t[1], t[0])
+	e2 := Sub(t[2], t[0])
+	return Cross(e1, e2)
+}
+
+func randomTriangle(rnd *rand.Rand) Triangle {
+	return Triangle{
+		randomVec(rnd),
+		randomVec(rnd),
+		randomVec(rnd),
+	}
+}