From 174f7cf7255e000a3ab2b25be3836b443bd6876e Mon Sep 17 00:00:00 2001 From: Paul Mach Date: Tue, 10 Nov 2015 15:21:23 -0800 Subject: [PATCH] Quadtree implementation --- bound.go | 9 + clustering/quadtree/quadtree.go | 46 ---- clustering/quadtree/quadtree_test.go | 11 - define.go | 4 +- point.go | 5 + pointers.go | 12 + quadtree/README.md | 43 +++ quadtree/benchmarks_test.go | 183 +++++++++++++ quadtree/examples_test.go | 42 +++ quadtree/quadtree.go | 395 +++++++++++++++++++++++++++ quadtree/quadtree_test.go | 110 ++++++++ 11 files changed, 800 insertions(+), 60 deletions(-) delete mode 100644 clustering/quadtree/quadtree.go delete mode 100644 clustering/quadtree/quadtree_test.go create mode 100644 pointers.go create mode 100644 quadtree/README.md create mode 100644 quadtree/benchmarks_test.go create mode 100644 quadtree/examples_test.go create mode 100644 quadtree/quadtree.go create mode 100644 quadtree/quadtree_test.go diff --git a/bound.go b/bound.go index 9efc5c5..966b53b 100644 --- a/bound.go +++ b/bound.go @@ -166,6 +166,15 @@ func geoHashInt2ranges(hash int64, bits int) (float64, float64, float64, float64 return lngMin, lngMax, latMin, latMax } +// Set allows for the modification of the bound values in place. +func (b *Bound) Set(west, east, south, north float64) { + b.sw[0] = west + b.sw[1] = south + + b.ne[0] = east + b.ne[1] = north +} + // Extend grows the bound to include the new point. func (b *Bound) Extend(point *Point) *Bound { diff --git a/clustering/quadtree/quadtree.go b/clustering/quadtree/quadtree.go deleted file mode 100644 index 9e5d8c9..0000000 --- a/clustering/quadtree/quadtree.go +++ /dev/null @@ -1,46 +0,0 @@ -package quadtree - -import ( - "github.com/paulmach/go.geo" -) - -type Quadtree struct { - TreeExtent *Extent -} - -type QuadTreeNode struct { - Leaf bool - Nodes []*QuadTreeNode - Point *geo.Point -} - -type Extent struct { - TopLeft *geo.Point - BottomRight *geo.Point -} - -type visitFunction func(*geo.Point, float64, float64, float64, float64) bool - -func NewFromPoints(points *geo.PointSet) *Quadtree { - return &Quadtree{} -} - -func New() *Quadtree { - return &Quadtree{} -} - -func (q *Quadtree) Add(p *geo.Point) *Quadtree { - return nil -} - -func (q *Quadtree) Visit(visit visitFunction) { - return nil -} - -func (q *Quadtree) Find(p *geo.Point) *geo.Point { - return nil -} - -func (q *Quadtree) SetExtent(e *Extent) { - q.TreeExtent = Extent -} diff --git a/clustering/quadtree/quadtree_test.go b/clustering/quadtree/quadtree_test.go deleted file mode 100644 index b8a60fa..0000000 --- a/clustering/quadtree/quadtree_test.go +++ /dev/null @@ -1,11 +0,0 @@ -package quadtree - -import ( - "testing" -) - -func TestSingleNodeNoBounds(t *testing.T) { -} - -func TestLocateSinglePoint(t *testing.T) { -} diff --git a/define.go b/define.go index f19948e..87979aa 100644 --- a/define.go +++ b/define.go @@ -1,8 +1,6 @@ package geo -import ( - "math" -) +import "math" var epsilon = 1e-6 diff --git a/point.go b/point.go index 2f61f42..744db88 100644 --- a/point.go +++ b/point.go @@ -64,6 +64,11 @@ func NewPointFromGeoHashInt64(hash int64, bits int) *Point { return NewPoint((west+east)/2.0, (north+south)/2.0) } +// Point, so point implements the pointer interface on itself. +func (p *Point) Point() *Point { + return p +} + // Transform applies a given projection or inverse projection to the current point. func (p *Point) Transform(projector Projector) *Point { projector(p) diff --git a/pointers.go b/pointers.go new file mode 100644 index 0000000..2f3877c --- /dev/null +++ b/pointers.go @@ -0,0 +1,12 @@ +package geo + +// A Pointer is the interface for something that has a point. +type Pointer interface { + // Point should return the "center" or other canonical point + // for the object. The caller is expected to Clone + // the point if changes need to be make. + Point() *Point +} + +// TODO: add some functionality around sets of pointers, +// ie. `type PointerSlice []Pointer` diff --git a/quadtree/README.md b/quadtree/README.md new file mode 100644 index 0000000..73e7a41 --- /dev/null +++ b/quadtree/README.md @@ -0,0 +1,43 @@ +go.geo/quadtree +=============== + +Package quadtree implements a quadtree using rectangular partitions. +Each point exists in a unique node; if multiple points are in the same position, +some points may be stored on internal nodes rather than leaf nodes. +This implementation is based heavily off of the +[d3 implementation](https://github.com/mbostock/d3/wiki/Quadtree-Geom). + +## Examples + + func ExampleQuadtreeFind() { + r := rand.New(rand.NewSource(42)) // to make things reproducible + + qt := quadtree.New(geo.NewBound(0, 1, 0, 1)) + + // insert 1000 random points + for i := 0; i < 1000; i++ { + qt.Insert(geo.NewPoint(r.Float64(), r.Float64())) + } + + nearest := qt.Find(geo.NewPoint(0.5, 0.5)) + fmt.Printf("nearest: %+v\n", nearest) + + // Output: + // nearest: POINT(0.4930591659434973 0.5196585530161364) + } + + func ExampleQuadtreeInBound() { + r := rand.New(rand.NewSource(52)) // to make things reproducible + + qt := quadtree.New(geo.NewBound(0, 1, 0, 1)) + + // insert 1000 random points + for i := 0; i < 1000; i++ { + qt.Insert(geo.NewPoint(r.Float64(), r.Float64())) + } + + bounded := qt.InBound(geo.NewBound(0.5, 0.5, 0.5, 0.5).Pad(0.05)) + fmt.Printf("in bound: %v\n", len(bounded)) + // Output: + // in bound: 10 + } diff --git a/quadtree/benchmarks_test.go b/quadtree/benchmarks_test.go new file mode 100644 index 0000000..3db116f --- /dev/null +++ b/quadtree/benchmarks_test.go @@ -0,0 +1,183 @@ +package quadtree + +import ( + "math" + "math/rand" + "testing" + + "github.com/paulmach/go.geo" +) + +func BenchmarkInsert(b *testing.B) { + r := rand.New(rand.NewSource(22)) + qt := New(geo.NewBound(0, 1, 0, 1)) + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + qt.Insert(geo.NewPoint(r.Float64(), r.Float64())) + } +} + +func BenchmarkFromPointer50(b *testing.B) { + r := rand.New(rand.NewSource(32)) + pointers := make([]geo.Pointer, 0, 50) + for i := 0; i < 50; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + NewFromPointers(pointers) + } +} + +func BenchmarkFromPointer100(b *testing.B) { + r := rand.New(rand.NewSource(42)) + pointers := make([]geo.Pointer, 0, 100) + for i := 0; i < 100; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + NewFromPointers(pointers) + } +} + +func BenchmarkFromPointer500(b *testing.B) { + r := rand.New(rand.NewSource(52)) + pointers := make([]geo.Pointer, 0, 500) + for i := 0; i < 500; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + NewFromPointers(pointers) + } +} + +func BenchmarkFromPointer1000(b *testing.B) { + r := rand.New(rand.NewSource(62)) + pointers := make([]geo.Pointer, 0, 1000) + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + NewFromPointers(pointers) + } +} + +func BenchmarkRandomFind1000(b *testing.B) { + r := rand.New(rand.NewSource(42)) + + var pointers []geo.Pointer + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + qt := NewFromPointers(pointers) + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + qt.Find(geo.NewPoint(r.Float64(), r.Float64())) + } +} + +func BenchmarkRandomFind1000Naive(b *testing.B) { + r := rand.New(rand.NewSource(42)) + + var pointers []geo.Pointer + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + looking := geo.NewPoint(r.Float64(), r.Float64()) + + min := math.MaxFloat64 + var best geo.Pointer + for _, p := range pointers { + if d := looking.SquaredDistanceFrom(p.Point()); d < min { + min = d + best = p + } + } + + _ = best + } +} + +func BenchmarkRandomInBound1000(b *testing.B) { + r := rand.New(rand.NewSource(43)) + + var pointers []geo.Pointer + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + qt := NewFromPointers(pointers) + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + p := geo.NewPoint(r.Float64(), r.Float64()) + qt.InBound(geo.NewBoundFromPoints(p, p).Pad(0.1)) + } +} + +func BenchmarkRandomInBound1000Naive(b *testing.B) { + r := rand.New(rand.NewSource(43)) + + var pointers []geo.Pointer + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + b.ReportAllocs() + b.ResetTimer() + + var near []geo.Pointer + for i := 0; i < b.N; i++ { + p := geo.NewPoint(r.Float64(), r.Float64()) + b := geo.NewBoundFromPoints(p, p).Pad(0.1) + + near = near[:0] + for _, p := range pointers { + if b.Contains(p.Point()) { + near = append(near, p) + } + } + + _ = len(near) + } +} + +func BenchmarkRandomInBound1000Buf(b *testing.B) { + r := rand.New(rand.NewSource(43)) + + var pointers []geo.Pointer + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + + qt := NewFromPointers(pointers) + + var buf []geo.Pointer + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + p := geo.NewPoint(r.Float64(), r.Float64()) + buf = qt.InBound(geo.NewBoundFromPoints(p, p).Pad(0.1), buf) + } +} diff --git a/quadtree/examples_test.go b/quadtree/examples_test.go new file mode 100644 index 0000000..109f451 --- /dev/null +++ b/quadtree/examples_test.go @@ -0,0 +1,42 @@ +package quadtree_test + +import ( + "fmt" + "math/rand" + + "github.com/paulmach/go.geo" + "github.com/paulmach/go.geo/quadtree" +) + +func ExampleQuadtreeFind() { + r := rand.New(rand.NewSource(42)) // to make things reproducable + + qt := quadtree.New(geo.NewBound(0, 1, 0, 1)) + + // insert 1000 random points + for i := 0; i < 1000; i++ { + qt.Insert(geo.NewPoint(r.Float64(), r.Float64())) + } + + nearest := qt.Find(geo.NewPoint(0.5, 0.5)) + fmt.Printf("nearest: %+v\n", nearest) + + // Output: + // nearest: POINT(0.4930591659434973 0.5196585530161364) +} + +func ExampleQuadtreeInBound() { + r := rand.New(rand.NewSource(52)) // to make things reproducable + + qt := quadtree.New(geo.NewBound(0, 1, 0, 1)) + + // insert 1000 random points + for i := 0; i < 1000; i++ { + qt.Insert(geo.NewPoint(r.Float64(), r.Float64())) + } + + bounded := qt.InBound(geo.NewBound(0.5, 0.5, 0.5, 0.5).Pad(0.05)) + fmt.Printf("in bound: %v\n", len(bounded)) + // Output: + // in bound: 10 +} diff --git a/quadtree/quadtree.go b/quadtree/quadtree.go new file mode 100644 index 0000000..f933feb --- /dev/null +++ b/quadtree/quadtree.go @@ -0,0 +1,395 @@ +// Package quadtree implements a quadtree using rectangular partitions. +// Each point exists in a unique node; if multiple points are in the same position, +// some points may be stored on internal nodes rather than leaf nodes. +// This implementation is based heavily off of the d3 implementation: +// https://github.com/mbostock/d3/wiki/Quadtree-Geom +package quadtree + +import ( + "errors" + "math" + + "github.com/paulmach/go.geo" +) + +var ( + // ErrPointOutsideOfBounds is returned when trying to add a point + // to a quad tree and the point is outside the bounds used to create the tree. + ErrPointOutsideOfBounds = errors.New("quadtree: point outside of bounds") +) + +// Quadtree implements a two-dimensional recursive spatial subdivision +// of geo.Pointers. This implementation uses rectangular partitions. +type Quadtree struct { + // Threshold indicates the limit of how deep the quadtree can go. + // Points closer than this will essentially be put in the same leaf node and stop recursion. + // The correct value depends on the use case. The default is computed + // off the bounds to keep the tree at most 12 levels deep. So points that + // are closer than 1/4096 * max(bound.width, bound.height) will essentially be + // stored in the same leaf node. For optimal tree performance you want this to happen + // sometimes but not very often. + Threshold float64 + + bound *geo.Bound + root *node + freeNodes []node + freeIndex int +} + +// node represents a node of the quad tree. Each node stores a Value +// and has links to its 4 children +type node struct { + children [4]*node + internal bool + pointer geo.Pointer +} + +// New creates a new quadtree for the given bound. Added points +// must be within this bound. +func New(bound *geo.Bound) *Quadtree { + return &Quadtree{ + Threshold: math.Max(bound.Width(), bound.Height()) / float64(1<<12), + bound: bound, + } +} + +// NewFromPointSet creates a quadtree from a pointset. +// Copies the points into the quad tree. Modifying the points later +// will invalidate the quad tree and lead to unexpected result. +func NewFromPointSet(set *geo.PointSet) *Quadtree { + q := New(set.Bound()) + q.freeNodes = make([]node, set.Length(), set.Length()) + + ps := []geo.Point(*set) + for i := range ps { + q.Insert(&ps[i]) + } + + return q +} + +// NewFromPointers creates a quadtree from a set of pointers. +func NewFromPointers(points []geo.Pointer) *Quadtree { + if len(points) == 0 { + // This is kind of meaningless but is what will happen + // if using an empty pointset above. + return New(geo.NewBound(0, 0, 0, 0)) + } + + b := geo.NewBoundFromPoints(points[0].Point(), points[0].Point()) + for _, p := range points { + b.Extend(p.Point()) + } + + q := New(b) + q.freeNodes = make([]node, len(points), len(points)) + + for _, p := range points { + q.Insert(p) + } + + return q +} + +// Bound returns the bounds used for the quad tree. +func (q *Quadtree) Bound() *geo.Bound { + return q.bound +} + +// Insert puts an object into the quad tree, must be within the quadtree bounds. +// If the pointer returns nil, the point will be ignored. +// This function is not thread-safe, ie. multiple goroutines cannot insert into +// a single quadtree. +func (q *Quadtree) Insert(p geo.Pointer) error { + if p == nil { + return nil + } + + point := p.Point() + if point == nil { + return nil + } + + if !q.bound.Contains(point) { + return ErrPointOutsideOfBounds + } + + if q.root == nil { + q.root = &node{} + } + + q.insert(q.root, p, + q.bound.Left(), q.bound.Right(), + q.bound.Bottom(), q.bound.Top(), + ) + + return nil +} + +// nextNode returns the next node from a preallocated list. +// This resulted in about 15% improvement in quadtree creation. +func (q *Quadtree) nextNode() *node { + if l := len(q.freeNodes); q.freeIndex >= l { + // Exponentially decrease the preallocation size. + // On a handful of tests, number of nodes was about 1.5 times pointers. + l /= 2 + + // min size of the preallocation. I think this could be bigger as it's + // not that much memory overhead. Optimizing this more would need + // to be use case specific. + if l < 25 { + l = 25 + } + + q.freeNodes = make([]node, l, l) + q.freeIndex = 0 + } + + n := &q.freeNodes[q.freeIndex] + q.freeIndex++ + return n +} + +func (q *Quadtree) insert(n *node, p geo.Pointer, left, right, bottom, top float64) { + point := p.Point() + if n.internal { + i := 0 + + // figure which child of this internal node the point is in. + if cy := (bottom + top) / 2.0; point.Y() <= cy { + top = cy + i = 2 + } else { + bottom = cy + } + + if cx := (left + right) / 2.0; point.X() >= cx { + left = cx + i++ + } else { + right = cx + } + + if n.children[i] == nil { + // child not yet created so automatically add the pointer to it and return. + n.children[i] = q.nextNode() + n.children[i].pointer = p + return + } + + // proceed down to the child to see if it's a leaf yet and we can add the pointer there. + q.insert(n.children[i], p, left, right, bottom, top) + return + } + + if n.pointer == nil { + // leaf without a pointer. I believe this only happens for the first pointer added. + // ie. initialized empty root node with no data. + n.pointer = p + return + } + + // leaf node with a point in it. Now we're splitting it and making it an internal node. + nPoint := n.pointer.Point() + n.internal = true + + if dx, dy := nPoint.X()-point.X(), nPoint.Y()-point.Y(); dx*dx+dy*dy <= q.Threshold*q.Threshold { + // similar/duplicate point to stop recursion. + i := childIndex((left+right)/2.0, (bottom+top)/2.0, point) + if n.children[i] == nil { + n.children[i] = q.nextNode() + n.children[i].pointer = p + return + } + q.insert(n, p, left, right, bottom, top) + return + } + + nPointer := n.pointer + n.pointer = nil + + // current node is now an internal node. + // first re add its pointer as one of its children, + // then add the new pointer as one of the children. + q.insert(n, nPointer, left, right, bottom, top) + q.insert(n, p, left, right, bottom, top) +} + +// Find returns the closest Value/Pointer in the quadtree. +// This function is thread safe. Multiple goroutines can read from +// a pre-created tree. +func (q *Quadtree) Find(p *geo.Point) geo.Pointer { + if q.root == nil { + return nil + } + + v := &findVisitor{ + point: p, + closestBound: q.bound.Clone(), + minDistSquared: math.MaxFloat64, + } + + newVisit(v).Visit(q.root, + q.bound.Left(), q.bound.Right(), + q.bound.Bottom(), q.bound.Top(), + ) + + return v.closest +} + +// InBound returns a slice with all the pointers in the quadtree that are +// within the given bound. An optional buffer parameter is provided to allow +// for the reuse of result slice memory. This function is thread safe. +// Multiple goroutines can read from a pre-created tree. +func (q *Quadtree) InBound(b *geo.Bound, buf ...[]geo.Pointer) []geo.Pointer { + if q.root == nil { + return nil + } + + var p []geo.Pointer + if len(buf) > 0 { + p = buf[0][:0] + } + v := &inBoundVisitor{ + bound: b, + pointers: p, + } + + newVisit(v).Visit(q.root, + q.bound.Left(), q.bound.Right(), + q.bound.Bottom(), q.bound.Top(), + ) + + return v.pointers +} + +// The visit stuff is a more go like (hopefully) implementation of the +// d3.quadtree.visit function. It is not exported, but if there is a +// good use case, it could be. + +type visitor interface { + // Bound returns the current relevant bound so we can prune irrelevant nodes + // from the search. + Bound() *geo.Bound + Visit(p geo.Pointer) + + // Point should return the specific point being search for, or null if there + // isn't one (ie. searching by bound). This helps guide the search to the + // best child node first. + Point() *geo.Point +} + +// visit provides a framework for walking the quad tree. +// Currently used by the `Find` and `InBound` functions. +type visit struct { + visitor visitor +} + +func newVisit(v visitor) *visit { + return &visit{ + visitor: v, + } +} + +func (v *visit) Visit(n *node, left, right, bottom, top float64) { + b := v.visitor.Bound() + if left > b.Right() || right < b.Left() || + bottom > b.Top() || top < b.Bottom() { + return + } + + if n.pointer != nil { + v.visitor.Visit(n.pointer) + } + + if !n.internal { + return + } + + cx := (left + right) / 2.0 + cy := (bottom + top) / 2.0 + + i := 0 + if p := v.visitor.Point(); p != nil { + // go to the child node the point is in first. + i = childIndex(cx, cy, p) + } + + for j := i; j < i+4; j++ { + if n.children[j%4] == nil { + continue + } + + if k := j % 4; k == 0 { + v.Visit(n.children[0], left, cx, cy, top) + } else if k == 1 { + v.Visit(n.children[1], cx, right, cy, top) + } else if k == 2 { + v.Visit(n.children[2], left, cx, bottom, cy) + } else if k == 3 { + v.Visit(n.children[3], cx, right, bottom, cy) + } + } +} + +type findVisitor struct { + point *geo.Point + closest geo.Pointer + closestBound *geo.Bound + minDistSquared float64 +} + +func (v *findVisitor) Bound() *geo.Bound { + return v.closestBound +} + +func (v *findVisitor) Point() *geo.Point { + return v.point +} + +func (v *findVisitor) Visit(p geo.Pointer) { + point := p.Point() + if d := point.SquaredDistanceFrom(v.point); d < v.minDistSquared { + v.minDistSquared = d + v.closest = p + + d = math.Sqrt(d) + x := v.point.X() + y := v.point.Y() + v.closestBound.Set(x-d, x+d, y-d, y+d) + } + + return +} + +type inBoundVisitor struct { + bound *geo.Bound + pointers []geo.Pointer +} + +func (v *inBoundVisitor) Bound() *geo.Bound { + return v.bound +} + +func (v *inBoundVisitor) Point() *geo.Point { + return nil +} + +func (v *inBoundVisitor) Visit(p geo.Pointer) { + if v.bound.Contains(p.Point()) { + v.pointers = append(v.pointers, p) + } +} + +func childIndex(cx, cy float64, point *geo.Point) int { + i := 0 + if point.Y() <= cy { + i = 2 + } + + if point.X() >= cx { + i++ + } + + return i +} diff --git a/quadtree/quadtree_test.go b/quadtree/quadtree_test.go new file mode 100644 index 0000000..e34f51a --- /dev/null +++ b/quadtree/quadtree_test.go @@ -0,0 +1,110 @@ +package quadtree + +import ( + "math/rand" + "testing" + + "github.com/paulmach/go.geo" +) + +func TestNew(t *testing.T) { + bound := geo.NewBound(0, 1, 2, 3) + qt := New(bound) + + if qt.Bound() != bound { + t.Errorf("should use provided bound, got %v", qt.Bound) + } + + ps := geo.NewPointSet() + ps.Push(geo.NewPoint(0, 2)) + ps.Push(geo.NewPoint(1, 3)) + + qt = NewFromPointSet(ps) + if !qt.Bound().Equals(ps.Bound()) { + t.Errorf("should take bound from pointset, got %v", qt.Bound()) + } +} + +func TestQuadtreeFind(t *testing.T) { + points := []geo.Pointer{} + dim := 17 + for i := 0; i < dim*dim; i++ { + points = append(points, geo.NewPoint(float64(i%dim), float64(i/dim))) + } + + q := NewFromPointers(points) + + // table test + type findTest struct { + Point *geo.Point + Expected *geo.Point + } + + tests := []findTest{ + {Point: geo.NewPoint(0.1, 0.1), Expected: geo.NewPoint(0, 0)}, + {Point: geo.NewPoint(3.1, 2.9), Expected: geo.NewPoint(3, 3)}, + {Point: geo.NewPoint(7.5, 7.5), Expected: geo.NewPoint(7, 7)}, + {Point: geo.NewPoint(7.1, 7.1), Expected: geo.NewPoint(7, 7)}, + {Point: geo.NewPoint(8.5, 8.5), Expected: geo.NewPoint(7, 7)}, + {Point: geo.NewPoint(0.1, 15.9), Expected: geo.NewPoint(0, 16)}, + {Point: geo.NewPoint(15.9, 15.9), Expected: geo.NewPoint(16, 16)}, + } + + for i, test := range tests { + if i != 3 { + continue + } + if v := q.Find(test.Point); !v.Point().Equals(test.Expected) { + t.Errorf("incorrect point on %d, got %v", i, v) + } + } +} + +func TestQuadtreeFindRandom(t *testing.T) { + r := rand.New(rand.NewSource(42)) + ps := geo.NewPointSet() + for i := 0; i < 1000; i++ { + ps.Push(geo.NewPoint(r.Float64(), r.Float64())) + } + qt := NewFromPointSet(ps) + + for i := 0; i < 1000; i++ { + p := geo.NewPoint(r.Float64(), r.Float64()) + + f := qt.Find(p) + _, j := ps.DistanceFrom(p) + + if e := ps.GetAt(j); !e.Equals(f.Point()) { + t.Errorf("index: %d, unexpected point %v != %v", i, e, f.Point()) + } + } +} + +func TestQuadtreeInBoundRandom(t *testing.T) { + r := rand.New(rand.NewSource(43)) + + var pointers []geo.Pointer + for i := 0; i < 1000; i++ { + pointers = append(pointers, geo.NewPoint(r.Float64(), r.Float64())) + } + qt := NewFromPointers(pointers) + + for i := 0; i < 1000; i++ { + p := geo.NewPoint(r.Float64(), r.Float64()) + + b := geo.NewBoundFromPoints(p, p).Pad(0.1) + ps := qt.InBound(b) + + // find the right answer brute force + var list []geo.Pointer + for _, p := range pointers { + if b.Contains(p.Point()) { + list = append(list, p) + } + } + + if len(list) != len(ps) { + t.Errorf("index: %d, lengths not equal %v != %v", i, len(list), len(ps)) + } + } +}