Mercurial > gemma

changeset 2488:cb55d7eaaa36 octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Started with the idea to clip an octree by an bounding polygon.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 03 Mar 2019 14:45:41 +0100
parents bd46ffbb944e
children 4f292ff74d9e
files cmd/octreediff/main.go pkg/octree/polygon.go pkg/octree/tree.go
diffstat 3 files changed, 178 insertions(+), 16 deletions(-) [+]
line wrap: on
line diff
--- a/cmd/octreediff/main.go	Sun Mar 03 13:00:50 2019 +0100
+++ b/cmd/octreediff/main.go	Sun Mar 03 14:45:41 2019 +0100
@@ -65,6 +65,23 @@
     4326
   )
 `
+	clippingPolygonSQL = `
+WITH joined AS (
+  SELECT
+    sr.area      AS area,
+    sr.date_info AS date_info
+  FROM waterway.sounding_results sr JOIN
+     waterway.bottlenecks bn ON sr.bottleneck_id = bn.id
+  WHERE bn.bottleneck_id = $1
+)
+SELECT ST_AsBinary(
+  ST_Buffer(
+  ST_intersection(
+    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
+    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
+  ),
+  0.001)) AS area
+`
 	insertContourSQLClipped = `
 WITH joined AS (
   SELECT
@@ -347,6 +364,21 @@
 		last = now
 		log.Printf("num points: %d\n", len(result))
 
+		var clip []byte
+
+		if err := tx.QueryRowContext(
+			ctx, clippingPolygonSQL,
+			bottleneck,
+			first.EPSG,
+			firstDate, secondDate,
+		).Scan(&clip); err != nil {
+			return err
+		}
+
+		now = time.Now()
+		log.Printf("loading clipping polygon took %v\n", now.Sub(last))
+		last = now
+
 		tri, err := result.triangulate()
 		if err != nil {
 			return err
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/polygon.go	Sun Mar 03 14:45:41 2019 +0100
@@ -0,0 +1,36 @@
+// This is Free Software under GNU Affero General Public License v >= 3.0
+// without warranty, see README.md and license for details.
+//
+// SPDX-License-Identifier: AGPL-3.0-or-later
+// License-Filename: LICENSES/AGPL-3.0.txt
+//
+// Copyright (C) 2018 by via donau
+//   – Österreichische Wasserstraßen-Gesellschaft mbH
+// Software engineering by Intevation GmbH
+//
+// Author(s):
+//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
+
+package octree
+
+type Polygon struct {
+	// TODO: Implement me!
+}
+
+type IntersectionType byte
+
+const (
+	IntersectionInside IntersectionType = iota
+	IntersectionOutSide
+	IntersectionOverlaps
+)
+
+func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType {
+	// TODO: Implement me
+	return IntersectionOutSide
+}
+
+func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType {
+	// TODO: Implement me
+	return IntersectionOutSide
+}
--- a/pkg/octree/tree.go	Sun Mar 03 13:00:50 2019 +0100
+++ b/pkg/octree/tree.go	Sun Mar 03 14:45:41 2019 +0100
@@ -32,6 +32,11 @@
 	Max Vertex
 }
 
+type boxFrame struct {
+	pos int32
+	Box2D
+}
+
 func (ot *Tree) Vertices() []Vertex {
 	return ot.vertices
 }
@@ -43,6 +48,105 @@
 	{0.5, 0.5, 1.0, 1.0},
 }
 
+func (ot *Tree) Clip(p *Polygon) {
+
+	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
+
+	stack := []boxFrame{{1, all}}
+
+	triChecks := make(map[int32]IntersectionType)
+
+frames:
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		if top.pos > 0 { // node
+			switch p.IntersectionBox2D(top.Box2D) {
+			case IntersectionInside:
+				// all inside so nothing to clip.
+				continue frames
+			case IntersectionOutSide:
+				// all outside -> clip from tree.
+				if index := ot.index[top.pos:]; len(index) > 7 {
+					for i := range index {
+						index[i] = 0
+					}
+				}
+				continue frames
+			default: // Overlaps
+				if index := ot.index[top.pos:]; len(index) > 7 {
+				children:
+					for i := 0; i < 4; i++ {
+						a := index[i]
+						b := index[i+4]
+						if a == 0 && b == 0 {
+							continue
+						}
+						dx := top.X2 - top.X1
+						dy := top.Y2 - top.Y1
+						nbox := Box2D{
+							dx*scale[i][0] + top.X1,
+							dy*scale[i][1] + top.Y1,
+							dx*scale[i][2] + top.X1,
+							dy*scale[i][3] + top.Y1,
+						}
+						switch p.IntersectionBox2D(nbox) {
+						case IntersectionInside:
+							// all inside so nothing to clip.
+							continue children
+						case IntersectionOutSide:
+							// all are ouside -> clip from tree.
+							index[i] = 0
+							index[i+4] = 0
+							continue
+						default: // Overlaps
+							if a != 0 {
+								stack = append(stack, boxFrame{a, nbox})
+							}
+							if b != 0 {
+								stack = append(stack, boxFrame{b, nbox})
+							}
+						}
+					}
+				}
+			}
+		} else { // leaf
+			pos := -top.pos - 1
+			n := ot.index[pos]
+			indices := ot.index[pos+1 : pos+1+n]
+		tris:
+			for i := len(indices) - 1; i >= 0; i-- {
+				triIndex := indices[i]
+				what, found := triChecks[triIndex]
+				if !found {
+					tri := ot.triangles[triIndex]
+					t := Triangle{
+						ot.vertices[tri[0]],
+						ot.vertices[tri[1]],
+						ot.vertices[tri[2]],
+					}
+					what = p.IntersectionWithTriangle(&t)
+					triChecks[triIndex] = what
+				}
+				switch what {
+				case IntersectionInside:
+					// triangle inside -> stay.
+					continue tris
+				default:
+					// outside or not fully covered -> remove.
+					if i < len(indices)-1 {
+						copy(indices[i:], indices[i+1:])
+					}
+					indices[len(indices)-1] = 0
+					indices = indices[:len(indices)-1]
+				}
+			}
+			ot.index[pos] = int32(len(indices))
+		}
+	}
+}
+
 func (ot *Tree) Value(x, y float64) (float64, bool) {
 
 	// out of bounding box
@@ -51,14 +155,9 @@
 		return 0, false
 	}
 
-	type frame struct {
-		pos int32
-		Box2D
-	}
-
 	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
 
-	stack := []frame{{1, all}}
+	stack := []boxFrame{{1, all}}
 
 	for len(stack) > 0 {
 		top := stack[len(stack)-1]
@@ -82,10 +181,10 @@
 					}
 					if nbox.Contains(x, y) {
 						if a != 0 {
-							stack = append(stack, frame{a, nbox})
+							stack = append(stack, boxFrame{a, nbox})
 						}
 						if b != 0 {
-							stack = append(stack, frame{b, nbox})
+							stack = append(stack, boxFrame{b, nbox})
 						}
 						break
 					}
@@ -131,18 +230,13 @@
 
 	line := NewPlane2D(x1, y1, x2, y2)
 
-	type frame struct {
-		pos int32
-		Box2D
-	}
-
 	dupes := map[int32]struct{}{}
 
 	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
 	//log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1))
 	//log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1))
 
-	stack := []frame{{1, all}}
+	stack := []boxFrame{{1, all}}
 
 	for len(stack) > 0 {
 		top := stack[len(stack)-1]
@@ -166,10 +260,10 @@
 					}
 					if nbox.Intersects(box) && nbox.IntersectsPlane(line) {
 						if a != 0 {
-							stack = append(stack, frame{a, nbox})
+							stack = append(stack, boxFrame{a, nbox})
 						}
 						if b != 0 {
-							stack = append(stack, frame{b, nbox})
+							stack = append(stack, boxFrame{b, nbox})
 						}
 					}
 				}