changeset 2490:c9164ff98871 octree-diff

Started with point in polygon check.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 03 Mar 2019 21:10:07 +0100
parents 4f292ff74d9e
children 10681749371d
files pkg/octree/polygon.go
diffstat 1 files changed, 77 insertions(+), 9 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/octree/polygon.go	Sun Mar 03 18:08:40 2019 +0100
+++ b/pkg/octree/polygon.go	Sun Mar 03 21:10:07 2019 +0100
@@ -22,9 +22,11 @@
 	"gemma.intevation.de/gemma/pkg/wkb"
 )
 
+type ring []float64
+
 type Polygon struct {
 	// TODO: Implement me!
-	rings [][]float64
+	rings []ring
 }
 
 type IntersectionType byte
@@ -45,6 +47,71 @@
 	return IntersectionOutSide
 }
 
+func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 }
+
+func (rng ring) contains(point []float64) bool {
+	if !rng.isClosed() {
+		return false
+	}
+
+	contains := intersectsWithRaycast(point, rng[:2], rng[len(rng)-2:len(rng)])
+
+	for i := 2; i < len(rng); i += 2 {
+		if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) {
+			contains = !contains
+		}
+	}
+
+	return contains
+}
+
+// Using the raycast algorithm, this returns whether or not the passed in point
+// Intersects with the edge drawn by the passed in start and end points.
+// Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go
+func intersectsWithRaycast(point, start, end []float64) bool {
+
+	// Always ensure that the the first point
+	// has a y coordinate that is less than the second point
+	if start[1] > end[1] {
+		// Switch the points if otherwise.
+		start, end = end, start
+	}
+
+	// Move the point's y coordinate
+	// outside of the bounds of the testing region
+	// so we can start drawing a ray
+	for point[1] == start[1] || point[1] == end[1] {
+		y := math.Nextafter(point[1], math.Inf(1))
+		point = []float64{point[0], y}
+	}
+
+	// If we are outside of the polygon, indicate so.
+	if point[1] < start[1] || point[1] > end[1] {
+		return false
+	}
+
+	if start[0] > end[0] {
+		if point[0] > start[0] {
+			return false
+		}
+		if point[0] < end[0] {
+			return true
+		}
+	} else {
+		if point[0] > end[0] {
+			return false
+		}
+		if point[0] < start[0] {
+			return true
+		}
+	}
+
+	raySlope := (point[1] - start[1]) / (point[0] - start[0])
+	diagSlope := (end[1] - start[1]) / (end[0] - start[0])
+
+	return raySlope >= diagSlope
+}
+
 func (p *Polygon) FromWKB(data []byte) error {
 
 	r := bytes.NewReader(data)
@@ -79,17 +146,18 @@
 		return err
 	}
 
-	rings := make([][]float64, numRings)
+	rngs := make([]ring, numRings)
 
-	for ring := uint32(0); ring < numRings; ring++ {
+	for rng := uint32(0); rng < numRings; rng++ {
 		var numVertices uint32
 		if err = binary.Read(r, order, &numVertices); err != nil {
 			return err
 		}
 
-		vertices := make([]float64, 2*numVertices)
+		numVertices *= 2
+		vertices := make([]float64, numVertices)
 
-		for v := uint32(0); v < numVertices; v++ {
+		for v := uint32(0); v < numVertices; v += 2 {
 			var lat, lon uint64
 			if err = binary.Read(r, order, &lat); err != nil {
 				return err
@@ -97,14 +165,14 @@
 			if err = binary.Read(r, order, &lon); err != nil {
 				return err
 			}
-			vertices[v*2] = math.Float64frombits(lat)
-			vertices[v*2+1] = math.Float64frombits(lon)
+			vertices[v] = math.Float64frombits(lat)
+			vertices[v+1] = math.Float64frombits(lon)
 		}
 
-		rings[ring] = vertices
+		rngs[rng] = vertices
 	}
 
-	p.rings = rings
+	p.rings = rngs
 
 	return nil
 }