Mercurial > gemma

changeset 4737:e8df3aeb0a32

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged stack-polygons branch back into default.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 18 Oct 2019 11:20:33 +0200
parents 7a40a39853a9 (current diff) fcbd585ae308 (diff)
children a0a4298a4756
files
diffstat 5 files changed, 278 insertions(+), 67 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/controllers/diff.go	Fri Oct 18 10:45:25 2019 +0200
+++ b/pkg/controllers/diff.go	Fri Oct 18 11:20:33 2019 +0200
@@ -84,15 +84,8 @@
   $2,
   ST_Transform(
     ST_Multi(
-      ST_CollectionExtract(
-        ST_MakeValid(
-          ST_Multi(
-             ST_Collectionextract(
-                ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5),
-                3
-             )
-            )
-        ),
+      ST_Collectionextract(
+        ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5),
         3
       )
     ),
--- a/pkg/imports/sr.go	Fri Oct 18 10:45:25 2019 +0200
+++ b/pkg/imports/sr.go	Fri Oct 18 11:20:33 2019 +0200
@@ -174,15 +174,8 @@
   $2,
   ST_Transform(
     ST_Multi(
-      ST_CollectionExtract(
-        ST_MakeValid(
-          ST_Multi(
-             ST_Collectionextract(
-                ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5),
-                3
-             )
-            )
-        ),
+      ST_Collectionextract(
+        ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5),
         3
       )
     ),
--- a/pkg/octree/areas.go	Fri Oct 18 10:45:25 2019 +0200
+++ b/pkg/octree/areas.go	Fri Oct 18 11:20:33 2019 +0200
@@ -189,6 +189,10 @@
 	reprojX := common.Linear(0.5, min.X, 1.5, min.X+cellSize)
 	reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize)
 
+	reproj := func(x, y float64) (float64, float64) {
+		return reprojX(x), reprojY(y)
+	}
+
 	cnts := make(chan int)
 
 	doContours := func() {
@@ -200,28 +204,7 @@
 				continue
 			}
 
-			// We need to bring it back to the
-			// none raster coordinate system.
-			a := make(wkb.MultiPolygonGeom, len(c))
-
-			for i, pl := range c {
-				shell := make(wkb.LinearRingGeom, len(pl))
-				for j, pt := range pl {
-					shell[j] = wkb.PointGeom{
-						X: reprojX(pt.X),
-						Y: reprojY(pt.Y),
-					}
-				}
-				/*
-					if !shell.CCW() {
-						log.Println("not ccw")
-						shell.Reverse()
-					}
-				*/
-				a[i] = wkb.PolygonGeom{shell}
-			}
-
-			areas[hIdx] = a
+			areas[hIdx] = buildMultipolygon(c, reproj)
 		}
 	}
 
@@ -240,3 +223,220 @@
 
 	return areas
 }
+
+type contour []contourmap.Point
+
+type bboxNode struct {
+	box      Box2D
+	cnt      contour
+	children []*bboxNode
+}
+
+type bboxForest struct {
+	roots []*bboxNode
+}
+
+func (cnt contour) contains(o contour) bool {
+	return contains(cnt, o[0].X, o[0].Y) ||
+		contains(cnt, o[len(o)/2].X, o[len(o)/2].Y)
+}
+
+func (cnt contour) closed() bool {
+	return len(cnt) >= 3
+}
+
+func (cnt contour) length() int {
+	return len(cnt)
+}
+
+func (cnt contour) point(i int) (float64, float64) {
+	return cnt[i].X, cnt[i].Y
+}
+
+func (cnt contour) bbox() Box2D {
+	minX, minY := math.MaxFloat64, math.MaxFloat64
+	maxX, maxY := -math.MaxFloat64, -math.MaxFloat64
+
+	for _, p := range cnt {
+		if p.X < minX {
+			minX = p.X
+		}
+		if p.X > maxX {
+			maxX = p.X
+		}
+		if p.Y < minY {
+			minY = p.Y
+		}
+		if p.Y > maxY {
+			maxY = p.Y
+		}
+	}
+	return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY}
+}
+
+func (bn *bboxNode) insert(cnt contour, box Box2D) {
+	// check if children are inside new
+	var nr *bboxNode
+
+	for i, r := range bn.children {
+		if r.box.Inside(box) && cnt.contains(r.cnt) {
+			if nr == nil {
+				nr = &bboxNode{box: box, cnt: cnt}
+			}
+			nr.children = append(nr.children, r)
+			bn.children[i] = nil
+		}
+	}
+
+	// we have a new child
+	if nr != nil {
+		// compact the list
+		for i := len(bn.children) - 1; i >= 0; i-- {
+			if bn.children[i] == nil {
+				if i < len(bn.children)-1 {
+					copy(bn.children[i:], bn.children[i+1:])
+				}
+				bn.children[len(bn.children)-1] = nil
+				bn.children = bn.children[:len(bn.children)-1]
+			}
+		}
+		bn.children = append(bn.children, nr)
+		return
+	}
+
+	// check if new is inside an old
+	for _, r := range bn.children {
+		if box.Inside(r.box) && r.cnt.contains(cnt) {
+			r.insert(cnt, box)
+			return
+		}
+	}
+
+	// its a new child node.
+	nr = &bboxNode{box: box, cnt: cnt}
+	bn.children = append(bn.children, nr)
+}
+
+func (bf *bboxForest) insert(cnt contour) {
+	box := cnt.bbox()
+
+	// check if roots are inside new
+	var nr *bboxNode
+
+	for i, r := range bf.roots {
+		if r.box.Inside(box) && cnt.contains(r.cnt) {
+			if nr == nil {
+				nr = &bboxNode{box: box, cnt: cnt}
+			}
+			nr.children = append(nr.children, r)
+			bf.roots[i] = nil
+		}
+	}
+
+	// we have a new root
+	if nr != nil {
+		// compact the list
+		for i := len(bf.roots) - 1; i >= 0; i-- {
+			if bf.roots[i] == nil {
+				if i < len(bf.roots)-1 {
+					copy(bf.roots[i:], bf.roots[i+1:])
+				}
+				bf.roots[len(bf.roots)-1] = nil
+				bf.roots = bf.roots[:len(bf.roots)-1]
+			}
+		}
+		bf.roots = append(bf.roots, nr)
+		return
+	}
+
+	// check if new is inside an old
+	for _, r := range bf.roots {
+		if box.Inside(r.box) && r.cnt.contains(cnt) {
+			r.insert(cnt, box)
+			return
+		}
+	}
+
+	// its a new root node.
+	nr = &bboxNode{box: box, cnt: cnt}
+	bf.roots = append(bf.roots, nr)
+}
+
+type bboxOutFunc func(contour, []contour)
+
+func (bn *bboxNode) generate(out bboxOutFunc) {
+	if len(bn.children) == 0 {
+		out(bn.cnt, nil)
+		return
+	}
+
+	var grands []*bboxNode
+
+	holes := make([]contour, len(bn.children))
+
+	for i, ch := range bn.children {
+		holes[i] = ch.cnt
+		grands = append(grands, ch.children...)
+	}
+	out(bn.cnt, holes)
+	for _, grand := range grands {
+		grand.generate(out)
+	}
+}
+
+func (bf *bboxForest) generate(out bboxOutFunc) {
+	for _, r := range bf.roots {
+		r.generate(out)
+	}
+}
+
+func buildMultipolygon(
+	cnts []contourmap.Contour,
+	reproj func(float64, float64) (float64, float64),
+) wkb.MultiPolygonGeom {
+
+	var bf bboxForest
+
+	for _, cnt := range cnts {
+		bf.insert(contour(cnt))
+	}
+
+	//log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots))
+
+	var mp wkb.MultiPolygonGeom
+
+	out := func(sh contour, hls []contour) {
+
+		polygon := make(wkb.PolygonGeom, 1+len(hls))
+
+		// Handle shell
+		shell := make(wkb.LinearRingGeom, len(sh))
+		for j, pt := range sh {
+			x, y := reproj(pt.X, pt.Y)
+			shell[j] = wkb.PointGeom{X: x, Y: y}
+		}
+		if shell.CCW() {
+			shell.Reverse()
+		}
+		polygon[0] = shell
+
+		// handle holes
+		for i, hl := range hls {
+			hole := make(wkb.LinearRingGeom, len(hl))
+			for j, pt := range hl {
+				x, y := reproj(pt.X, pt.Y)
+				hole[j] = wkb.PointGeom{X: x, Y: y}
+			}
+			if !hole.CCW() {
+				hole.Reverse()
+			}
+			polygon[1+i] = hole
+		}
+
+		mp = append(mp, polygon)
+	}
+
+	bf.generate(out)
+
+	return mp
+}
--- a/pkg/octree/polygon.go	Fri Oct 18 10:45:25 2019 +0200
+++ b/pkg/octree/polygon.go	Fri Oct 18 11:20:33 2019 +0200
@@ -281,19 +281,18 @@
 
 	// No intersection -> check inside or outside
 	// if an abritrary point  is inside or not.
-	point := []float64{box.X1, box.Y1}
 
 	// Check holes first: inside a hole means outside.
 	if len(p.rings) > 1 {
 		for _, hole := range p.rings[1:] {
-			if hole.contains(point) {
+			if contains(hole, box.X1, box.Y1) {
 				return IntersectionOutSide
 			}
 		}
 	}
 
 	// Check shell
-	if p.rings[0].contains(point) {
+	if contains(p.rings[0], box.X1, box.Y1) {
 		return IntersectionInside
 	}
 	return IntersectionOutSide
@@ -324,37 +323,59 @@
 	}
 	// No intersection -> check inside or outside
 	// if an abritrary point  is inside or not.
-	point := []float64{t[0].X, t[0].Y}
+	pX, pY := t[0].X, t[0].Y
 
 	// Check holes first: inside a hole means outside.
 	if len(p.rings) > 1 {
 		for _, hole := range p.rings[1:] {
-			if hole.contains(point) {
+			if contains(hole, pX, pY) {
 				return IntersectionOutSide
 			}
 		}
 	}
 
 	// Check shell
-	if p.rings[0].contains(point) {
+	if contains(p.rings[0], pX, pY) {
 		return IntersectionInside
 	}
 	return IntersectionOutSide
 }
 
-func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 }
+func (rng ring) closed() bool {
+	return (len(rng) / 2) >= 3
+}
+
+func (rng ring) length() int {
+	return len(rng) / 2
+}
 
-func (rng ring) contains(point []float64) bool {
-	if !rng.isClosed() {
+func (rng ring) point(i int) (float64, float64) {
+	i *= 2
+	return rng[i], rng[i+1]
+}
+
+type segments interface {
+	closed() bool
+	length() int
+	point(int) (float64, float64)
+}
+
+func contains(s segments, pX, pY float64) bool {
+	if !s.closed() {
 		return false
 	}
 
-	end := len(rng)/2 - 1
+	n := s.length()
+
+	sX, sY := s.point(0)
+	eX, eY := s.point(n - 1)
 
-	contains := intersectsWithRaycast(point, rng[:2], rng[end*2:end*2+2])
+	contains := intersectsWithRaycast(pX, pY, sX, sY, eX, eY)
 
-	for i := 2; i < len(rng); i += 2 {
-		if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) {
+	for i := 1; i < n; i++ {
+		sX, sY := s.point(i - 1)
+		eX, eY := s.point(i)
+		if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) {
 			contains = !contains
 		}
 	}
@@ -365,46 +386,45 @@
 // 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 {
+func intersectsWithRaycast(pX, pY, sX, sY, eX, eY 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] {
+	if sY > eY {
 		// Switch the points if otherwise.
-		start, end = end, start
+		sX, sY, eX, eY = eX, eY, sX, sY
 	}
 
 	// 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}
+	for pY == sY || pY == eY {
+		pY = math.Nextafter(pY, math.Inf(1))
 	}
 
 	// If we are outside of the polygon, indicate so.
-	if point[1] < start[1] || point[1] > end[1] {
+	if pY < sY || pY > eY {
 		return false
 	}
 
-	if start[0] > end[0] {
-		if point[0] > start[0] {
+	if sX > eX {
+		if pX > sX {
 			return false
 		}
-		if point[0] < end[0] {
+		if pX < eX {
 			return true
 		}
 	} else {
-		if point[0] > end[0] {
+		if pX > eX {
 			return false
 		}
-		if point[0] < start[0] {
+		if pX < sX {
 			return true
 		}
 	}
 
-	raySlope := (point[1] - start[1]) / (point[0] - start[0])
-	diagSlope := (end[1] - start[1]) / (end[0] - start[0])
+	raySlope := (pY - sY) / (pX - sX)
+	diagSlope := (eY - sY) / (eX - sX)
 
 	return raySlope >= diagSlope
 }
--- a/pkg/octree/vertex.go	Fri Oct 18 10:45:25 2019 +0200
+++ b/pkg/octree/vertex.go	Fri Oct 18 11:20:33 2019 +0200
@@ -108,6 +108,11 @@
 	}
 }
 
+func (a Box2D) Inside(b Box2D) bool {
+	return a.X1 >= b.X1 && a.X2 <= b.X2 &&
+		a.Y1 >= b.Y1 && a.Y2 <= b.Y2
+}
+
 func (p Plane3D) Z(x, y float64) float64 {
 	// p.A*x + p.B*y + p.C*z + p.D = 0
 	return -(p.A*x + p.B*y + p.D) / p.C