spatial/r3: add Box

spatial/r3/box.go

@@ -0,0 +1,120 @@
+// 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"
+
+// Box is a 3D bounding box. Well formed Boxes Min components
+// are smaller than Max components.
+type Box struct {
+	Min, Max Vec
+}
+
+// NewBox is shorthand for Box{Min:Vec{x0,y0,z0}, Max:Vec{x1,y1,z1}}.
+// The sides are swapped so that the resulting Box is well formed.
+func NewBox(x0, y0, z0, x1, y1, z1 float64) Box {
+	return Box{
+		Min: Vec{X: math.Min(x0, x1), Y: math.Min(y0, y1), Z: math.Min(z0, z1)},
+		Max: Vec{X: math.Max(x0, x1), Y: math.Max(y0, y1), Z: math.Max(z0, z1)},
+	}
+}
+
+// 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 || a.Min.Z >= a.Max.Z
+}
+
+// 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))
+}
+
+// Vertices returns a slice of the 8 vertices
+// corresponding to each of the Box's corners.
+//
+// Ordering of vertices 0-3 is CCW in the XY plane starting at box minimum.
+// Ordering of vertices 4-7 is CCW in the XY plane starting at box minimum
+// for X and Y values and maximum Z value.
+//
+// Edges for the box can be constructed with the following indices:
+//  edges := [12][2]int{
+//   {0, 1}, {1, 2}, {2, 3}, {3, 0},
+//   {4, 5}, {5, 6}, {6, 7}, {7, 4},
+//   {0, 4}, {1, 5}, {2, 6}, {3, 7},
+//  }
+func (a Box) Vertices() []Vec {
+	return []Vec{
+		0: a.Min,
+		1: {X: a.Max.X, Y: a.Min.Y, Z: a.Min.Z},
+		2: {X: a.Max.X, Y: a.Max.Y, Z: a.Min.Z},
+		3: {X: a.Min.X, Y: a.Max.Y, Z: a.Min.Z},
+		4: {X: a.Min.X, Y: a.Min.Y, Z: a.Max.Z},
+		5: {X: a.Max.X, Y: a.Min.Y, Z: a.Max.Z},
+		6: a.Max,
+		7: {X: a.Min.X, Y: a.Max.Y, Z: a.Max.Z},
+	}
+}
+
+// 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 dimension
+// is scaled by scale.X. Negative elements of scale are interpreted as zero.
+func (a Box) Scale(scale Vec) Box {
+	scale = maxElem(scale, Vec{})
+	// TODO(soypat): Probably a better way to do this.
+	return centeredBox(a.Center(), mulElem(scale, a.Size()))
+}
+
+// centeredBox creates a Box with a given center and size.
+// Negative components of size will be interpreted as zero.
+func centeredBox(center, size Vec) Box {
+	size = maxElem(size, Vec{}) // set negative values to zero.
+	half := Scale(0.5, 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 &&
+		a.Min.Z <= v.Z && v.Z <= a.Max.Z
+}
+
+// Canon returns the canonical version of a. The returned Box has minimum
+// and maximum coordinates swapped if necessary so that it is well-formed.
+func (a Box) Canon() Box {
+	return Box{
+		Min: minElem(a.Min, a.Max),
+		Max: maxElem(a.Min, a.Max),
+	}
+}

spatial/r3/box_test.go

@@ -0,0 +1,200 @@
+// 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 (
+	"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 := [6]Vec{
+			Add(b.Max, Vec{1, 0, 0}),
+			Add(b.Max, Vec{0, 1, 0}),
+			Add(b.Max, Vec{0, 0, 1}),
+			Sub(b.Min, Vec{1, 0, 0}),
+			Sub(b.Min, Vec{0, 1, 0}),
+			Sub(b.Min, Vec{0, 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 := [6]Vec{
+			Add(maxElem(b1.Max, b2.Max), Vec{1, 0, 0}),
+			Add(maxElem(b1.Max, b2.Max), Vec{0, 1, 0}),
+			Add(maxElem(b1.Max, b2.Max), Vec{0, 0, 1}),
+			Sub(minElem(b1.Min, b2.Min), Vec{1, 0, 0}),
+			Sub(minElem(b1.Min, b2.Min), Vec{0, 1, 0}),
+			Sub(minElem(b1.Min, b2.Min), Vec{0, 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 !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 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 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 8, got %d", len(gots))
+		}
+		for j, want := range wants {
+			got := gots[j]
+			if !vecEqual(want, got) {
+				t.Errorf("%dth vertex not equal", j)
+			}
+		}
+	}
+}
+
+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")
+		}
+	}
+}
+
+// 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),
+		Z: randomRange(b.Min.Z, b.Max.Z),
+	}
+}
+
+// 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: {X: a.Max.X, Y: a.Min.Y, Z: a.Min.Z},
+		2: {X: a.Max.X, Y: a.Max.Y, Z: a.Min.Z},
+		3: {X: a.Min.X, Y: a.Max.Y, Z: a.Min.Z},
+		4: {X: a.Min.X, Y: a.Min.Y, Z: a.Max.Z},
+		5: {X: a.Max.X, Y: a.Min.Y, Z: a.Max.Z},
+		6: a.Max,
+		7: {X: a.Min.X, Y: a.Max.Y, Z: a.Max.Z},
+	}
+}

spatial/r3/vector.go

@@ -83,7 +83,50 @@ func Cos(p, q Vec) float64 {
 	return Dot(p, q) / (Norm(p) * Norm(q))
 }
 
-// Box is a 3D bounding box.
+// minElem return a vector with the minimum components of two vectors.
-type Box struct {
+func minElem(a, b Vec) Vec {
-	Min, Max Vec
+	return Vec{
+		X: math.Min(a.X, b.X),
+		Y: math.Min(a.Y, b.Y),
+		Z: math.Min(a.Z, b.Z),
+	}
+}
+
+// maxElem return a vector with the maximum components of two vectors.
+func maxElem(a, b Vec) Vec {
+	return Vec{
+		X: math.Max(a.X, b.X),
+		Y: math.Max(a.Y, b.Y),
+		Z: math.Max(a.Z, b.Z),
+	}
+}
+
+// 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),
+		Z: math.Abs(a.Z),
+	}
+}
+
+// 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,
+		Z: a.Z * b.Z,
+	}
+}
+
+// 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,
+		Z: a.Z / b.Z,
+	}
 }