diff --git a/form2/must2/basic.go b/form2/must2/basic.go
index 8a258a1..1e1abcc 100644
--- a/form2/must2/basic.go
+++ b/form2/must2/basic.go
@@ -31,7 +31,7 @@ func (s *circle) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a 2d circle.
-func (s *circle) BoundingBox() r2.Box {
+func (s *circle) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -60,7 +60,7 @@ func (s *box) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a 2d box.
-func (s *box) BoundingBox() r2.Box {
+func (s *box) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -92,6 +92,6 @@ func (s *line) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a 2d line.
-func (s *line) BoundingBox() r2.Box {
+func (s *line) Bounds() r2.Box {
 	return s.bb
 }
diff --git a/form2/must2/poly.go b/form2/must2/poly.go
index 45d5bf7..bc80389 100644
--- a/form2/must2/poly.go
+++ b/form2/must2/poly.go
@@ -107,7 +107,7 @@ func (s *polygon) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a 2d polygon.
-func (s *polygon) BoundingBox() r2.Box {
+func (s *polygon) Bounds() r2.Box {
 	return s.bb
 }
 
diff --git a/form3/basic.go b/form3/basic.go
index a14a7b0..ed37617 100644
--- a/form3/basic.go
+++ b/form3/basic.go
@@ -80,8 +80,8 @@ func Cone(height, r0, r1, round float64) (s sdf.SDF3, err error) {
 // ChamferedCylinder intersects a chamfered cylinder with an SDF3.
 func ChamferedCylinder(s sdf.SDF3, kb, kt float64) (sdf.SDF3, error) {
 	// get the length and radius from the bounding box
-	l := s.BoundingBox().Max.Z
-	r := s.BoundingBox().Max.X
+	l := s.Bounds().Max.Z
+	r := s.Bounds().Max.X
 	p := form2.NewPolygon()
 	p.Add(0, -l)
 	p.Add(r, -l).Chamfer(r * kb)
diff --git a/form3/must3/cylinders.go b/form3/must3/cylinders.go
index ed11dba..1f5bd29 100644
--- a/form3/must3/cylinders.go
+++ b/form3/must3/cylinders.go
@@ -41,7 +41,7 @@ func (s *box) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a 3d box.
-func (s *box) BoundingBox() r3.Box {
+func (s *box) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -72,7 +72,7 @@ func (s *sphere) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a sphere.
-func (s *sphere) BoundingBox() r3.Box {
+func (s *sphere) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -121,7 +121,7 @@ func (s *cylinder) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a cylinder.
-func (s *cylinder) BoundingBox() r3.Box {
+func (s *cylinder) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -201,7 +201,7 @@ func (s *cone) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox return the bounding box for the trucated cone..
-func (s *cone) BoundingBox() r3.Box {
+func (s *cone) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -234,8 +234,8 @@ func sdfBox3d(p, s r3.Vec) float64 {
 // ChamferedCylinder intersects a chamfered cylinder with an SDF3.
 func ChamferedCylinder(s sdf.SDF3, kb, kt float64) sdf.SDF3 {
 	// get the length and radius from the bounding box
-	l := s.BoundingBox().Max.Z
-	r := s.BoundingBox().Max.X
+	l := s.Bounds().Max.Z
+	r := s.Bounds().Max.X
 	p := form2.NewPolygon()
 	p.Add(0, -l)
 	p.Add(r, -l).Chamfer(r * kb)
diff --git a/form3/must3/screw.go b/form3/must3/screw.go
index 6aa6340..661d1b8 100644
--- a/form3/must3/screw.go
+++ b/form3/must3/screw.go
@@ -79,6 +79,6 @@ func (s *screw) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a 3d screw form.
-func (s *screw) BoundingBox() r3.Box {
+func (s *screw) Bounds() r3.Box {
 	return s.bb
 }
diff --git a/matrix.go b/matrix.go
index 20deb56..c3bbc4a 100644
--- a/matrix.go
+++ b/matrix.go
@@ -510,3 +510,47 @@ func (a m22) Inverse() m22 {
 	m.x11 = a.x00 * d
 	return m
 }
+
+// rotateToVector returns the rotation matrix that transforms a onto the same direction as b.
+func rotateToVec(a, b r3.Vec) m44 {
+	// is either vector == 0?
+	if d3.EqualWithin(a, r3.Vec{}, epsilon) || d3.EqualWithin(b, r3.Vec{}, epsilon) {
+		return Identity3d()
+	}
+	// normalize both vectors
+	a = r3.Unit(a)
+	b = r3.Unit(b)
+	// are the vectors the same?
+	if d3.EqualWithin(a, b, epsilon) {
+		return Identity3d()
+	}
+
+	// are the vectors opposite (180 degrees apart)?
+	if d3.EqualWithin(r3.Scale(-1, a), b, epsilon) {
+		return m44{
+			-1, 0, 0, 0,
+			0, -1, 0, 0,
+			0, 0, -1, 0,
+			0, 0, 0, 1,
+		}
+	}
+	// general case
+	// See:	https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
+	v := r3.Cross(a, b)
+	vx := r3.Skew(v)
+
+	k := 1 / (1 + r3.Dot(a, b))
+	vx2 := r3.NewMat(nil)
+	vx2.Mul(vx, vx)
+	vx2.Scale(k, vx2)
+
+	// Calculate sum of matrices.
+	vx.Add(vx, r3.Eye())
+	vx.Add(vx, vx2)
+	return m44{
+		vx.At(0, 0), vx.At(0, 1), vx.At(0, 2), 0,
+		vx.At(1, 0), vx.At(1, 1), vx.At(1, 2), 0,
+		vx.At(2, 0), vx.At(2, 1), vx.At(2, 2), 0,
+		0, 0, 0, 1,
+	}
+}
diff --git a/render/internal_test.go b/render/internal_test.go
index 1230390..3200cc8 100644
--- a/render/internal_test.go
+++ b/render/internal_test.go
@@ -35,7 +35,7 @@ func TestSTLWriteReadback(t *testing.T) {
 		TotalLength: 40.,
 		ShankLength: 10.0,
 	})
-	size := r3.Norm(d3.Box(s0.BoundingBox()).Size())
+	size := r3.Norm(d3.Box(s0.Bounds()).Size())
 	// calculate relative tolerance
 	rtol := tol * size / quality
 	input, err := RenderAll(NewOctreeRenderer(s0, quality))
diff --git a/render/kdrender.go b/render/kdrender.go
index 47dd655..fa7ef9b 100644
--- a/render/kdrender.go
+++ b/render/kdrender.go
@@ -72,7 +72,7 @@ func (s kdSDF) Nearest(v r3.Vec) kdTriangle {
 	return got.(kdTriangle)
 }
 
-func (s kdSDF) BoundingBox() r3.Box {
+func (s kdSDF) Bounds() r3.Box {
 	bb := s.tree.Root.Bounding
 	if bb == nil {
 		panic("got nil bounding box?")
diff --git a/render/kdrender_test.go b/render/kdrender_test.go
index 382449c..be781d2 100644
--- a/render/kdrender_test.go
+++ b/render/kdrender_test.go
@@ -29,7 +29,7 @@ func TestKDSDF(t *testing.T) {
 	}
 	t.Error(len(model), "triangles")
 	kdf := render.NewKDSDF(model)
-	t.Error(kdf.BoundingBox())
+	t.Error(kdf.Bounds())
 	start := time.Now()
 	outside := kdf.Evaluate(r3.Vec{2, 0, 0}) // evaluate point outside bounds
 	inside := kdf.Evaluate(r3.Vec{0, 0, 0})  // evaluate point inside bounds
diff --git a/render/octree_renderer.go b/render/octree_renderer.go
index 617e63d..b426d2e 100644
--- a/render/octree_renderer.go
+++ b/render/octree_renderer.go
@@ -35,7 +35,7 @@ func NewOctreeRenderer(s sdf.SDF3, meshCells int) *octree {
 	}
 	// Scale the bounding box about the center to make sure the boundaries
 	// aren't on the object surface.
-	bb := d3.Box(s.BoundingBox())
+	bb := d3.Box(s.Bounds())
 	bb = bb.ScaleAboutCenter(1.01)
 	longAxis := d3.Max(bb.Size())
 	// We want to test the smallest cube (side == resolution) for emptiness
diff --git a/sdf2.go b/sdf2.go
index 7627ebb..86c4f4f 100644
--- a/sdf2.go
+++ b/sdf2.go
@@ -1,9 +1,3 @@
-/*
-
-2D Signed Distance Functions
-
-*/
-
 package sdf
 
 import (
@@ -16,10 +10,17 @@ import (
 	"gonum.org/v1/gonum/spatial/r3"
 )
 
+// 2D signed distance function utility functions.
+
 // SDF2 is the interface to a 2d signed distance function object.
 type SDF2 interface {
+	// Evaluate takes a point in 2D space as input and returns
+	// the minimum distance of the SDF2 to the point. The distance
+	// is negative if the point is contained within the SDF3.
 	Evaluate(p r2.Vec) float64
-	BoundingBox() r2.Box
+
+	// Bounds returns the bounding box that completely contains the SDF2.
+	Bounds() r2.Box
 }
 
 type SDF2Union interface {
@@ -69,7 +70,7 @@ func Cut2D(sdf SDF2, a, v r2.Vec) SDF2 {
 	v = r2.Unit(v)
 	s.n = r2.Vec{-v.Y, v.X}
 	// TODO - cut the bounding box
-	s.bb = sdf.BoundingBox()
+	s.bb = sdf.Bounds()
 	return &s
 }
 
@@ -79,7 +80,7 @@ func (s *CutSDF2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for the cut SDF2.
-func (s *CutSDF2) BoundingBox() r2.Box {
+func (s *CutSDF2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -98,7 +99,7 @@ func Transform2D(sdf SDF2, m m33) SDF2 {
 	s := TransformSDF2{}
 	s.sdf = sdf
 	s.mInv = m.Inverse()
-	s.bb = m.MulBox(sdf.BoundingBox())
+	s.bb = m.MulBox(sdf.Bounds())
 	return &s
 }
 
@@ -110,7 +111,7 @@ func (s *TransformSDF2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a transformed SDF2.
-func (s *TransformSDF2) BoundingBox() r2.Box {
+func (s *TransformSDF2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -131,7 +132,7 @@ func ScaleUniform2D(sdf SDF2, k float64) SDF2 {
 		sdf:  sdf,
 		k:    k,
 		invk: 1.0 / k,
-		bb:   m.MulBox(sdf.BoundingBox()),
+		bb:   m.MulBox(sdf.Bounds()),
 	}
 }
 
@@ -142,20 +143,20 @@ func (s *ScaleUniformSDF2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of an SDF2 with uniform scaling.
-func (s *ScaleUniformSDF2) BoundingBox() r2.Box {
+func (s *ScaleUniformSDF2) Bounds() r2.Box {
 	return s.bb
 }
 
 // Center2D centers the origin of an SDF2 on it's bounding box.
 func Center2D(s SDF2) SDF2 {
-	ofs := r2.Scale(-1, d2.Box(s.BoundingBox()).Center())
+	ofs := r2.Scale(-1, d2.Box(s.Bounds()).Center())
 	return Transform2D(s, Translate2d(ofs))
 }
 
 // CenterAndScale2D centers the origin of an SDF2 on it's bounding box, and then scales it.
 // Distance is correct with scaling.
 func CenterAndScale2D(s SDF2, k float64) SDF2 {
-	ofs := r2.Scale(-1, d2.Box(s.BoundingBox()).Center())
+	ofs := r2.Scale(-1, d2.Box(s.Bounds()).Center())
 	s = Transform2D(s, Translate2d(ofs))
 	return ScaleUniform2D(s, k)
 }
@@ -183,7 +184,7 @@ func Array2D(sdf SDF2, num V2i, step r2.Vec) SDF2Union {
 	s.step = step
 	s.min = math.Min
 	// work out the bounding box
-	bb0 := d2.Box(sdf.BoundingBox())
+	bb0 := d2.Box(sdf.Bounds())
 	// TODO verify
 	bb1 := bb0.Translate(d2.MulElem(step, r2.Sub(R2FromI(num), d2.Elem(1)))) // step.Mul(num.SubScalar(1).Tor2.Vec()))
 	s.bb = r2.Box(bb0.Extend(bb1))
@@ -208,7 +209,7 @@ func (s *array2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a grid array of SDF2s.
-func (s *array2) BoundingBox() r2.Box {
+func (s *array2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -233,7 +234,7 @@ func RotateUnion2D(sdf SDF2, num int, step m33) SDF2 {
 	s.step = step.Inverse()
 	s.min = math.Min
 	// work out the bounding box
-	vset := d2.Box(sdf.BoundingBox()).Vertices()
+	vset := d2.Box(sdf.Bounds()).Vertices()
 	bbMin := vset[0]
 	bbMax := vset[0]
 	for i := 0; i < s.num; i++ {
@@ -264,7 +265,7 @@ func (s *rotateUnion2) SetMin(min MinFunc) {
 }
 
 // BoundingBox returns the bounding box of a union of rotated SDF2s.
-func (s *rotateUnion2) BoundingBox() r2.Box {
+func (s *rotateUnion2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -285,7 +286,7 @@ func RotateCopy2D(sdf SDF2, n int) SDF2 {
 	s.sdf = sdf
 	s.theta = 2 * math.Pi / float64(n)
 	// work out the bounding box
-	bb := d2.Box(sdf.BoundingBox())
+	bb := d2.Box(sdf.Bounds())
 	rmax := 0.0
 	// find the bounding box vertex with the greatest distance from the origin
 	for _, v := range bb.Vertices() {
@@ -306,7 +307,7 @@ func (s *rotateCopy2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a rotate/copy SDF2.
-func (s *rotateCopy2) BoundingBox() r2.Box {
+func (s *rotateCopy2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -342,7 +343,7 @@ func Slice2D(sdf SDF3, a, n r3.Vec) SDF2 {
 	// TODO: This is bigger than it needs to be. We could consider intersection
 	// between the plane and the edges of the 3d bounding box for a smaller 2d
 	// bounding box in some circumstances.
-	v3 := d3.Box(sdf.BoundingBox()).Vertices()
+	v3 := d3.Box(sdf.Bounds()).Vertices()
 	vec := make(d2.Set, len(v3))
 	n = r3.Unit(n)
 	for i, v := range v3 {
@@ -364,7 +365,7 @@ func (s *slice2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of the sliced SDF2.
-func (s *slice2) BoundingBox() r2.Box {
+func (s *slice2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -387,9 +388,9 @@ func Union2D(sdf ...SDF2) SDF2Union {
 		}
 	}
 	// work out the bounding box
-	bb := d2.Box(s.sdf[0].BoundingBox())
+	bb := d2.Box(s.sdf[0].Bounds())
 	for _, x := range s.sdf {
-		bb = bb.Extend(d2.Box(x.BoundingBox()))
+		bb = bb.Extend(d2.Box(x.Bounds()))
 	}
 	s.bb = r2.Box(bb)
 	s.min = math.Min
@@ -403,7 +404,7 @@ func (s *union2) Evaluate(p r2.Vec) float64 {
 	minDist2 := -1.0
 	minIndex := 0
 	for i := range s.sdf {
-		vs[i] = d2.Box(s.sdf[i].BoundingBox()).MinMaxDist2(p)
+		vs[i] = d2.Box(s.sdf[i].Bounds()).MinMaxDist2(p)
 		// as we go record the sdf with the minimum minimum d2 value
 		if minDist2 < 0 || vs[i].X < minDist2 {
 			minDist2 = vs[i].X
@@ -449,7 +450,7 @@ func (s *union2) SetMin(min MinFunc) {
 }
 
 // BoundingBox returns the bounding box of an SDF2 union.
-func (s *union2) BoundingBox() r2.Box {
+func (s *union2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -470,7 +471,7 @@ func Difference2D(s0, s1 SDF2) SDF2Diff {
 	s.s0 = s0
 	s.s1 = s1
 	s.max = math.Max
-	s.bb = s0.BoundingBox()
+	s.bb = s0.Bounds()
 	return &s
 }
 
@@ -485,7 +486,7 @@ func (s *diff2) SetMax(max MaxFunc) {
 }
 
 // BoundingBox returns the bounding box of the difference of two SDF2s.
-func (s *diff2) BoundingBox() r2.Box {
+func (s *diff2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -505,7 +506,7 @@ func Elongate2D(sdf SDF2, h r2.Vec) SDF2 {
 		hn:  r2.Scale(0.5, h),
 	}
 	// bounding box
-	bb := d2.Box(sdf.BoundingBox())
+	bb := d2.Box(sdf.Bounds())
 	bb0 := bb.Translate(s.hp)
 	bb1 := bb.Translate(s.hn)
 	s.bb = r2.Box(bb0.Extend(bb1))
@@ -519,7 +520,7 @@ func (s *elongate2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of an elongated SDF2.
-func (s *elongate2) BoundingBox() r2.Box {
+func (s *elongate2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -527,7 +528,7 @@ func (s *elongate2) BoundingBox() r2.Box {
 func generateMesh2D(s SDF2, grid V2i) (d2.Set, error) {
 
 	// create the grid mapping for the bounding box
-	m, err := newMap2(d2.Box(s.BoundingBox()), grid, false)
+	m, err := newMap2(d2.Box(s.Bounds()), grid, false)
 	if err != nil {
 		return nil, err
 	}
@@ -596,7 +597,7 @@ func Offset2D(sdf SDF2, offset float64) SDF2 {
 	s.sdf = sdf
 	s.offset = offset
 	// work out the bounding box
-	bb := d2.Box(sdf.BoundingBox())
+	bb := d2.Box(sdf.Bounds())
 	s.bb = r2.Box(d2.NewBox2(bb.Center(), r2.Add(bb.Size(), d2.Elem(2*offset)))) //NewBox2(bb.Center(), r2.Add(bb.Size(), d2.Elem(2*offset)))
 	return &s
 }
@@ -607,7 +608,7 @@ func (s *offset2) Evaluate(p r2.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of an offset SDF2.
-func (s *offset2) BoundingBox() r2.Box {
+func (s *offset2) Bounds() r2.Box {
 	return s.bb
 }
 
@@ -629,7 +630,7 @@ func Intersect2D(s0, s1 SDF2) SDF2Diff {
 	s.s1 = s1
 	s.max = math.Max
 	// TODO fix bounding box
-	s.bb = s0.BoundingBox()
+	s.bb = s0.Bounds()
 	return &s
 }
 
@@ -644,12 +645,12 @@ func (s *intersection2) SetMax(max MaxFunc) {
 }
 
 // BoundingBox returns the bounding box of an SDF2 intersection.
-func (s *intersection2) BoundingBox() r2.Box {
+func (s *intersection2) Bounds() r2.Box {
 	return s.bb
 }
 func empty2From(s SDF2) empty2 {
 	return empty2{
-		center: d2.Box(s.BoundingBox()).Center(),
+		center: d2.Box(s.Bounds()).Center(),
 	}
 }
 
@@ -663,7 +664,7 @@ func (e empty2) Evaluate(r2.Vec) float64 {
 	return math.MaxFloat64
 }
 
-func (e empty2) BoundingBox() r2.Box {
+func (e empty2) Bounds() r2.Box {
 	return r2.Box{
 		Min: e.center,
 		Max: e.center,
diff --git a/sdf3.go b/sdf3.go
index e419e5d..85410fd 100644
--- a/sdf3.go
+++ b/sdf3.go
@@ -9,10 +9,17 @@ import (
 	"gonum.org/v1/gonum/spatial/r3"
 )
 
+// 3D signed distance utility functions.
+
 // SDF3 is the interface to a 3d signed distance function object.
 type SDF3 interface {
+	// Evaluate takes a point in 3D space as input and returns
+	// the minimum distance of the SDF3 to the point. The distance
+	// is negative if the point is contained within the SDF3.
 	Evaluate(p r3.Vec) float64
-	BoundingBox() r3.Box
+	// Bounds returns the bounding box that completely contains
+	// the SDF3.
+	Bounds() r3.Box
 }
 
 type SDF3Union interface {
@@ -97,7 +104,7 @@ func (s *revolution3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a solid of revolution.
-func (s *revolution3) BoundingBox() r3.Box {
+func (s *revolution3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -177,7 +184,7 @@ func (s *extrude3) SetExtrude(extrude ExtrudeFunc) {
 }
 
 // BoundingBox returns the bounding box for an extrusion.
-func (s *extrude3) BoundingBox() r3.Box {
+func (s *extrude3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -250,7 +257,7 @@ func (s *extrudeRounded) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a rounded extrusion.
-func (s *extrudeRounded) BoundingBox() r3.Box {
+func (s *extrudeRounded) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -327,7 +334,7 @@ func (s *loft3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box for a loft extrusion.
-func (s *loft3) BoundingBox() r3.Box {
+func (s *loft3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -347,7 +354,7 @@ func Transform3D(sdf SDF3, matrix m44) SDF3 {
 	s.sdf = sdf
 	s.matrix = matrix
 	s.inverse = matrix.Inverse()
-	s.bb = matrix.MulBox(sdf.BoundingBox())
+	s.bb = matrix.MulBox(sdf.Bounds())
 	return &s
 }
 
@@ -358,7 +365,7 @@ func (s *transform3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a transformed SDF3.
-func (s *transform3) BoundingBox() r3.Box {
+func (s *transform3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -378,7 +385,7 @@ func ScaleUniform3D(sdf SDF3, k float64) SDF3 {
 		sdf:  sdf,
 		k:    k,
 		invK: 1.0 / k,
-		bb:   m.MulBox(sdf.BoundingBox()),
+		bb:   m.MulBox(sdf.Bounds()),
 	}
 }
 
@@ -390,7 +397,7 @@ func (s *scaleUniform3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a uniformly scaled SDF3.
-func (s *scaleUniform3) BoundingBox() r3.Box {
+func (s *scaleUniform3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -417,9 +424,9 @@ func Union3D(sdf ...SDF3) SDF3Union {
 		}
 	}
 	// work out the bounding box
-	bb := d3.Box(s.sdf[0].BoundingBox())
+	bb := d3.Box(s.sdf[0].Bounds())
 	for _, x := range s.sdf {
-		bb = bb.Extend(d3.Box(x.BoundingBox()))
+		bb = bb.Extend(d3.Box(x.Bounds()))
 	}
 	s.bb = r3.Box(bb)
 	s.min = math.Min
@@ -445,7 +452,7 @@ func (s *union3) SetMin(min MinFunc) {
 }
 
 // BoundingBox returns the bounding box of an SDF3 union.
-func (s *union3) BoundingBox() r3.Box {
+func (s *union3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -467,7 +474,7 @@ func Difference3D(s0, s1 SDF3) SDF3Diff {
 	s.s0 = s0
 	s.s1 = s1
 	s.max = math.Max
-	s.bb = s0.BoundingBox()
+	s.bb = s0.Bounds()
 	return &s
 }
 
@@ -482,7 +489,7 @@ func (s *diff3) SetMax(max MaxFunc) {
 }
 
 // BoundingBox returns the bounding box of the SDF3 difference.
-func (s *diff3) BoundingBox() r3.Box {
+func (s *diff3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -502,7 +509,7 @@ func Elongate3D(sdf SDF3, h r3.Vec) SDF3 {
 		hn:  r3.Scale(-0.5, h),
 	}
 	// bounding box
-	bb := d3.Box(sdf.BoundingBox())
+	bb := d3.Box(sdf.Bounds())
 	bb0 := bb.Translate(s.hp)
 	bb1 := bb.Translate(s.hn)
 	s.bb = r3.Box(bb0.Extend(bb1))
@@ -516,7 +523,7 @@ func (s *elongate3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of an elongated SDF3.
-func (s *elongate3) BoundingBox() r3.Box {
+func (s *elongate3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -539,7 +546,7 @@ func Intersect3D(s0, s1 SDF3) SDF3Diff {
 	s.s1 = s1
 	s.max = math.Max
 	// TODO fix bounding box
-	s.bb = s0.BoundingBox()
+	s.bb = s0.Bounds()
 	return &s
 }
 
@@ -554,7 +561,7 @@ func (s *intersection3) SetMax(max MaxFunc) {
 }
 
 // BoundingBox returns the bounding box of an SDF3 intersection.
-func (s *intersection3) BoundingBox() r3.Box {
+func (s *intersection3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -574,7 +581,7 @@ func Cut3D(sdf SDF3, a, n r3.Vec) SDF3 {
 	s.a = a
 	s.n = r3.Scale(-1, r3.Unit(n))
 	// TODO - cut the bounding box
-	s.bb = sdf.BoundingBox()
+	s.bb = sdf.Bounds()
 	return &s
 }
 
@@ -584,7 +591,7 @@ func (s *cut3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of the cut SDF3.
-func (s *cut3) BoundingBox() r3.Box {
+func (s *cut3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -609,7 +616,7 @@ func Array3D(sdf SDF3, num V3i, step r3.Vec) SDF3Union {
 	s.step = step
 	s.min = math.Min
 	// work out the bounding box
-	bb0 := d3.Box(sdf.BoundingBox())
+	bb0 := d3.Box(sdf.Bounds())
 	bb1 := bb0.Translate(d3.MulElem(step, num.SubScalar(1).ToV3()))
 	s.bb = r3.Box(bb0.Extend(bb1))
 	return &s
@@ -635,7 +642,7 @@ func (s *array3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of an XYZ SDF3 array.
-func (s *array3) BoundingBox() r3.Box {
+func (s *array3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -661,7 +668,7 @@ func RotateUnion3D(sdf SDF3, num int, step m44) SDF3Union {
 	s.step = step.Inverse()
 	s.min = math.Min
 	// work out the bounding box
-	v := d3.Box(sdf.BoundingBox()).Vertices()
+	v := d3.Box(sdf.Bounds()).Vertices()
 	bbMin := v[0]
 	bbMax := v[0]
 	for i := 0; i < s.num; i++ {
@@ -692,7 +699,7 @@ func (s *rotateUnion) SetMin(min MinFunc) {
 }
 
 // BoundingBox returns the bounding box of a rotate/union object.
-func (s *rotateUnion) BoundingBox() r3.Box {
+func (s *rotateUnion) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -714,7 +721,7 @@ func RotateCopy3D(sdf SDF3, num int) SDF3 {
 	s.sdf = sdf
 	s.theta = tau / float64(num)
 	// work out the bounding box
-	bb := d3.Box(sdf.BoundingBox())
+	bb := d3.Box(sdf.Bounds())
 	zmax := bb.Max.Z
 	zmin := bb.Min.Z
 	rmax := 0.0
@@ -739,7 +746,7 @@ func (s *rotateCopy3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a rotate/copy SDF3.
-func (s *rotateCopy3) BoundingBox() r3.Box {
+func (s *rotateCopy3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -780,8 +787,8 @@ func (s *ConnectedSDF3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a connected SDF3.
-func (s *ConnectedSDF3) BoundingBox() d3.Box {
-	return s.sdf.BoundingBox()
+func (s *ConnectedSDF3) Bounds() d3.Box {
+	return s.sdf.Bounds()
 }
 
 */
@@ -800,7 +807,7 @@ func Offset3D(sdf SDF3, offset float64) SDF3 {
 		distance: offset,
 	}
 	// bounding box
-	bb := d3.Box(sdf.BoundingBox())
+	bb := d3.Box(sdf.Bounds())
 	s.bb = r3.Box(d3.NewBox(bb.Center(), r3.Add(bb.Size(), d3.Elem(2*offset))))
 	return &s
 }
@@ -811,7 +818,7 @@ func (s *offset3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of an offset SDF3.
-func (s *offset3) BoundingBox() r3.Box {
+func (s *offset3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -827,7 +834,7 @@ func Shell3D(sdf SDF3, thickness float64) SDF3 {
 	if thickness <= 0 {
 		return empty3From(sdf)
 	}
-	bb := d3.Box(sdf.BoundingBox())
+	bb := d3.Box(sdf.Bounds())
 	return &shell3{
 		sdf:   sdf,
 		delta: 0.5 * thickness,
@@ -841,7 +848,7 @@ func (s *shell3) Evaluate(p r3.Vec) float64 {
 }
 
 // BoundingBox returns the bounding box of a shelled SDF3.
-func (s *shell3) BoundingBox() r3.Box {
+func (s *shell3) Bounds() r3.Box {
 	return s.bb
 }
 
@@ -894,7 +901,7 @@ func Orient3D(s SDF3, base r3.Vec, directions d3.Set) SDF3 {
 
 func empty3From(s SDF3) empty3 {
 	return empty3{
-		center: d3.Box(s.BoundingBox()).Center(),
+		center: d3.Box(s.Bounds()).Center(),
 	}
 }
 
@@ -908,7 +915,7 @@ func (e empty3) Evaluate(r3.Vec) float64 {
 	return math.MaxFloat64
 }
 
-func (e empty3) BoundingBox() r3.Box {
+func (e empty3) Bounds() r3.Box {
 	return r3.Box{
 		Min: e.center,
 		Max: e.center,
diff --git a/vec3.go b/vec3.go
deleted file mode 100644
index 9a5a2ef..0000000
--- a/vec3.go
+++ /dev/null
@@ -1,50 +0,0 @@
-package sdf
-
-import (
-	"github.com/soypat/sdf/internal/d3"
-	"gonum.org/v1/gonum/spatial/r3"
-)
-
-// rotateToVector returns the rotation matrix that transforms a onto the same direction as b.
-func rotateToVec(a, b r3.Vec) m44 {
-	// is either vector == 0?
-	if d3.EqualWithin(a, r3.Vec{}, epsilon) || d3.EqualWithin(b, r3.Vec{}, epsilon) {
-		return Identity3d()
-	}
-	// normalize both vectors
-	a = r3.Unit(a)
-	b = r3.Unit(b)
-	// are the vectors the same?
-	if d3.EqualWithin(a, b, epsilon) {
-		return Identity3d()
-	}
-
-	// are the vectors opposite (180 degrees apart)?
-	if d3.EqualWithin(r3.Scale(-1, a), b, epsilon) {
-		return m44{
-			-1, 0, 0, 0,
-			0, -1, 0, 0,
-			0, 0, -1, 0,
-			0, 0, 0, 1,
-		}
-	}
-	// general case
-	// See:	https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
-	v := r3.Cross(a, b)
-	vx := r3.Skew(v)
-
-	k := 1 / (1 + r3.Dot(a, b))
-	vx2 := r3.NewMat(nil)
-	vx2.Mul(vx, vx)
-	vx2.Scale(k, vx2)
-
-	// Calculate sum of matrices.
-	vx.Add(vx, r3.Eye())
-	vx.Add(vx, vx2)
-	return m44{
-		vx.At(0, 0), vx.At(0, 1), vx.At(0, 2), 0,
-		vx.At(1, 0), vx.At(1, 1), vx.At(1, 2), 0,
-		vx.At(2, 0), vx.At(2, 1), vx.At(2, 2), 0,
-		0, 0, 0, 1,
-	}
-}
diff --git a/veci.go b/veci.go
index e1251f6..1868fdf 100644
--- a/veci.go
+++ b/veci.go
@@ -1,13 +1,9 @@
-/*
-
-Integer 2D/3D Vectors
-
-*/
-
 package sdf
 
 import "gonum.org/v1/gonum/spatial/r3"
 
+// Integer vectors
+
 // V2i is a 2D integer vector.
 type V2i [2]int