First pack of modules: point definition and distance calculation functions.

Author: smira

Makefile

@@ -0,0 +1,7 @@
+all: test
+
+test:
+	go test -v ./... -gocheck.v
+
+bench:
+	go test -v -bench=. ./... -gocheck.b

README.md

@@ -1,2 +1,5 @@
-# go-point-clustering
-(Lat, lon) points fast clustering using DBScan algorithm 
+# Point Clutering
+
+(Lat, lon) points fast clustering using DBScan algorithm in Golang.
+
+

distance.go

@@ -0,0 +1,67 @@
+package cluster
+
+import (
+	"math"
+)
+
+const (
+	// Translating from degrees to radians
+	DegreeRad = math.Pi / 180.0
+	// Earth radius
+	EarthR = 6371.0
+)
+
+// DistanceSheprical is a spherical (optimized) distance between two points
+//
+// Result is distance in kilometers
+func DistanceSpherical(p1, p2 *Point) float64 {
+	v1 := (p1[1] - p2[1]) * DegreeRad
+	v1 = v1 * v1
+
+	v2 := (p1[0] - p2[0]) * DegreeRad * math.Cos((p1[1]+p2[1])/2.0*DegreeRad)
+	v2 = v2 * v2
+
+	return EarthR * math.Sqrt(v1+v2)
+}
+
+// FastSine caclulates sinus approximated to parabola
+//
+// Taken from: http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648
+func FastSine(x float64) float64 {
+	const (
+		B = 4 / math.Pi
+		C = -4 / (math.Pi * math.Pi)
+		P = 0.225
+	)
+
+	if x > math.Pi || x < -math.Pi {
+		panic("out of range")
+	}
+
+	y := B*x + C*x*math.Abs(x)
+	return P*(y*math.Abs(y)-y) + y
+}
+
+// FastCos calculates cosinus from sinus
+func FastCos(x float64) float64 {
+	x += math.Pi / 2.0
+	if x > math.Pi {
+		x -= 2 * math.Pi
+	}
+
+	return FastSine(x)
+}
+
+// Spherical distance with fast cosine and w/o sqrt and
+// normalization to Earth radius/radians
+//
+// To get real distance in km, take sqrt and multiply result by EarthR*DegreeRad
+//
+// In this library eps (distance) is adjusted so that we don't need
+// to do sqrt and multiplication
+func DistanceSphericalFast(p1, p2 *Point) float64 {
+	v1 := (p1[1] - p2[1])
+	v2 := (p1[0] - p2[0]) * FastCos((p1[1]+p2[1])/2.0*DegreeRad)
+
+	return v1*v1 + v2*v2
+}

distance_test.go

@@ -0,0 +1,53 @@
+package cluster
+
+import (
+	. "gopkg.in/check.v1"
+	"math"
+)
+
+type DistanceSuite struct {
+	p1, p2 Point
+}
+
+var _ = Suite(&DistanceSuite{})
+
+func (s *DistanceSuite) SetUpTest(c *C) {
+	s.p1 = Point{30.244759, 59.955982}
+	s.p2 = Point{30.24472, 59.955975}
+}
+
+func (s *DistanceSuite) TestFastCos(c *C) {
+	c.Check(FastCos(0), Equals, math.Cos(0))
+	c.Check(math.Abs(FastCos(0.1)-math.Cos(0.1)) < 0.001, Equals, true)
+	c.Check(math.Abs(FastCos(-0.1)-math.Cos(-0.1)) < 0.001, Equals, true)
+	c.Check(math.Abs(FastCos(1.0)-math.Cos(1.0)) < 0.001, Equals, true)
+}
+
+func (s *DistanceSuite) TestDistanceSpherical(c *C) {
+	c.Check(DistanceSpherical(&s.p1, &s.p2), Equals, 0.0023064907653812116)
+	c.Check(DistanceSpherical(&s.p2, &s.p1), Equals, 0.0023064907653812116)
+	c.Check(DistanceSpherical(&s.p1, &s.p1), Equals, 0.0)
+	c.Check(DistanceSpherical(&s.p2, &s.p2), Equals, 0.0)
+}
+
+func (s *DistanceSuite) TestDistanceSphericalFast(c *C) {
+	c.Check(DistanceSphericalFast(&s.p1, &s.p2), Equals, 4.3026720164084415e-10)
+	c.Check(DistanceSphericalFast(&s.p2, &s.p1), Equals, 4.3026720164084415e-10)
+	c.Check(DistanceSphericalFast(&s.p1, &s.p1), Equals, 0.0)
+	c.Check(DistanceSphericalFast(&s.p2, &s.p2), Equals, 0.0)
+
+	c.Check(math.Abs(math.Sqrt(DistanceSphericalFast(&s.p1, &s.p2))*DegreeRad*EarthR-
+		DistanceSpherical(&s.p1, &s.p2)) < 0.000001, Equals, true)
+}
+
+func (s *DistanceSuite) BenchmarkDistanceSpherical(c *C) {
+	for i := 0; i < c.N; i++ {
+		DistanceSpherical(&s.p1, &s.p2)
+	}
+}
+
+func (s *DistanceSuite) BenchmarkDistanceSphericalFast(c *C) {
+	for i := 0; i < c.N; i++ {
+		DistanceSphericalFast(&s.p1, &s.p2)
+	}
+}

point.go

@@ -0,0 +1,65 @@
+// Package cluster implements DBScan clustering on (lat, lon) using K-D Tree
+package cluster
+
+// Point is longitue, latittude
+type Point [2]float64
+
+// PointList is a slice of Points
+type PointList []Point
+
+// Cluster is a result of DBScan work
+type Cluster struct {
+	C      int
+	Points []int
+}
+
+// sqDist returns squared (w/o sqrt & normalization) distance between two points
+func (a *Point) sqDist(b *Point) float64 {
+	return DistanceSphericalFast(a, b)
+}
+
+// LessEq - a <= b
+func (a *Point) LessEq(b *Point) bool {
+	return a[0] <= b[0] && a[1] <= b[1]
+}
+
+// GreaterEq - a >= b
+func (a *Point) GreaterEq(b *Point) bool {
+	return a[0] >= b[0] && a[1] >= b[1]
+}
+
+// CentroidAndBounds calculates center and cluster bounds
+func (c *Cluster) CentroidAndBounds(points PointList) (center, min, max Point) {
+	if len(c.Points) == 0 {
+		panic("empty cluster")
+	}
+
+	min = Point{180.0, 90.0}
+	max = Point{-180.0, -90.0}
+
+	for _, i := range c.Points {
+		pt := points[i]
+
+		for j := range pt {
+			center[j] += pt[j]
+
+			if pt[j] < min[j] {
+				min[j] = pt[j]
+			}
+			if pt[j] > max[j] {
+				max[j] = pt[j]
+			}
+		}
+	}
+
+	for j := range center {
+		center[j] /= float64(len(c.Points))
+	}
+
+	return
+}
+
+// Inside checks if (innerMin, innerMax) rectangle is inside (outerMin, outMax) rectangle
+func Inside(innerMin, innerMax, outerMin, outerMax *Point) bool {
+	return innerMin.GreaterEq(outerMin) && innerMax.LessEq(outerMax)
+}

point_test.go

@@ -0,0 +1,30 @@
+package cluster
+
+import (
+	. "gopkg.in/check.v1"
+	"testing"
+)
+
+// Launch gocheck tests
+func Test(t *testing.T) {
+	TestingT(t)
+}
+
+type PointSuite struct {
+	points PointList
+}
+
+var _ = Suite(&PointSuite{})
+
+func (s *PointSuite) SetUpTest(c *C) {
+	s.points = PointList{Point{30.244759, 59.955982}, Point{30.24472, 59.955975}, Point{30.244358, 59.96698}}
+}
+
+func (s *PointSuite) TestCentroidAndBounds(c *C) {
+	c1 := Cluster{C: 0, Points: []int{0, 1, 2}}
+
+	center, min, max := c1.CentroidAndBounds(s.points)
+	c.Check(center, DeepEquals, Point{30.244612333333333, 59.95964566666667})
+	c.Check(min, DeepEquals, Point{30.244358, 59.955975})
+	c.Check(max, DeepEquals, Point{30.244759, 59.96698})
+}