From e34e6b933b2bc85fed92e07705c6c005e3a27f98 Mon Sep 17 00:00:00 2001
From: Dan Kortschak <dan@kortschak.io>
Date: Thu, 20 Jun 2019 09:38:02 +0930
Subject: [PATCH] spatial/{barneshut,vptree}: refuse points at infinity

---
 spatial/barneshut/barneshut2.go | 11 +++++++++++
 spatial/barneshut/barneshut3.go | 11 +++++++++++
 spatial/vptree/vptree.go        | 18 ++++++++++++++++--
 3 files changed, 38 insertions(+), 2 deletions(-)

diff --git a/spatial/barneshut/barneshut2.go b/spatial/barneshut/barneshut2.go
index e229b654..d7738a5c 100644
--- a/spatial/barneshut/barneshut2.go
+++ b/spatial/barneshut/barneshut2.go
@@ -47,12 +47,23 @@ type Plane struct {
 }
 
 // NewPlane returns a new Plane.
+//
+// Points in p must not be infinitely distant, otherwise NewPlane will panic.
 func NewPlane(p []Particle2) *Plane {
+	for _, l := range p {
+		if isInfR2(l.Coord2()) {
+			panic("barneshut: point at infinity")
+		}
+	}
 	q := Plane{Particles: p}
 	q.Reset()
 	return &q
 }
 
+func isInfR2(v r2.Vec) bool {
+	return math.IsInf(v.X, 0) || math.IsInf(v.Y, 0)
+}
+
 // Reset reconstructs the Barnes-Hut tree. Reset must be called if the
 // Particles field or elements of Particles have been altered, unless
 // ForceOn is called with theta=0 or no data structures have been
diff --git a/spatial/barneshut/barneshut3.go b/spatial/barneshut/barneshut3.go
index c19cf773..3d1a5bc9 100644
--- a/spatial/barneshut/barneshut3.go
+++ b/spatial/barneshut/barneshut3.go
@@ -47,12 +47,23 @@ type Volume struct {
 }
 
 // NewVolume returns a new Volume.
+//
+// Points in p must not be infinitely distant, otherwise NewVolume will panic.
 func NewVolume(p []Particle3) *Volume {
+	for _, l := range p {
+		if isInfR3(l.Coord3()) {
+			panic("barneshut: point at infinity")
+		}
+	}
 	q := Volume{Particles: p}
 	q.Reset()
 	return &q
 }
 
+func isInfR3(v r3.Vec) bool {
+	return math.IsInf(v.X, 0) || math.IsInf(v.Y, 0) || math.IsInf(v.Z, 0)
+}
+
 // Reset reconstructs the Barnes-Hut tree. Reset must be called if the
 // Particles field or elements of Particles have been altered, unless
 // ForceOn is called with theta=0 or no data structures have been
diff --git a/spatial/vptree/vptree.go b/spatial/vptree/vptree.go
index afd05721..458c223c 100644
--- a/spatial/vptree/vptree.go
+++ b/spatial/vptree/vptree.go
@@ -64,6 +64,8 @@ type Tree struct {
 // The order of elements in p will be altered after New returns. The src parameter
 // provides the source of randomness for vantage point selection. If src is nil
 // global rand package functions are used.
+//
+// Points in p must not be infinitely distant, otherwise New will panic.
 func New(p []Comparable, effort int, src rand.Source) *Tree {
 	var intn func(int) int
 	var shuf func(n int, swap func(i, j int))
@@ -118,13 +120,21 @@ func (b *builder) selectVantage(s []Comparable, effort int) Comparable {
 	choices := b.random(effort, s)
 	for _, p := range choices {
 		for i, q := range choices {
-			b.work[i] = p.Distance(q)
+			d := p.Distance(q)
+			if math.IsInf(d, 0) {
+				panic("vptree: point at infinity")
+			}
+			b.work[i] = d
 		}
 		variance := stat.Variance(b.work, nil)
 		if variance > bestVar {
 			best, bestVar = p, variance
 		}
 	}
+	if best == nil {
+		// This should never be reached.
+		panic("vptree: could not find vantage point")
+	}
 	return best
 }
 
@@ -139,7 +149,11 @@ func (b *builder) random(n int, s []Comparable) []Comparable {
 func (b *builder) partition(v Comparable, s []Comparable) (radius float64, closer, further []Comparable) {
 	b.work = b.work[:len(s)]
 	for i, p := range s {
-		b.work[i] = v.Distance(p)
+		d := v.Distance(p)
+		if math.IsInf(d, 0) {
+			panic("vptree: point at infinity")
+		}
+		b.work[i] = d
 	}
 	sort.Sort(byDist{dists: b.work, points: s})