From c79d628c46a17df396ec692b940458133411ae1e Mon Sep 17 00:00:00 2001 From: soypat Date: Sun, 8 May 2022 11:54:09 -0300 Subject: [PATCH] spatial/r3: add Box --- spatial/r3/box.go | 120 +++++++++++++++++++++++++ spatial/r3/box_test.go | 200 +++++++++++++++++++++++++++++++++++++++++ spatial/r3/vector.go | 49 +++++++++- 3 files changed, 366 insertions(+), 3 deletions(-) create mode 100644 spatial/r3/box.go create mode 100644 spatial/r3/box_test.go diff --git a/spatial/r3/box.go b/spatial/r3/box.go new file mode 100644 index 00000000..e9130d65 --- /dev/null +++ b/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), + } +} diff --git a/spatial/r3/box_test.go b/spatial/r3/box_test.go new file mode 100644 index 00000000..2424ae5c --- /dev/null +++ b/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}, + } +} diff --git a/spatial/r3/vector.go b/spatial/r3/vector.go index da34c8ab..bcd0f881 100644 --- a/spatial/r3/vector.go +++ b/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. -type Box struct { - Min, Max Vec +// minElem return a vector with the minimum components of two vectors. +func minElem(a, b Vec) 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, + } }