Merge pull request #30 from paulmach/ml/quadtree

A quadtree of pointers

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 {
 

define.go

@@ -1,8 +1,6 @@
 package geo
 
-import (
-	"math"
-)
+import "math"
 
 var epsilon = 1e-6
 

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)

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`

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
+	}

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)
+	}
+}

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
+}