Mercurial > gemma

changeset 4658:4bbfe3dd2ab5 stree-experiment

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Completed usage of STRTrees.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 14 Oct 2019 14:58:04 +0200
parents a2f16987911b
children 92a7551d626f
files pkg/controllers/cross.go pkg/controllers/diff.go pkg/imports/isr.go pkg/octree/cache.go pkg/octree/loader.go pkg/octree/strtree.go pkg/octree/tree.go
diffstat 7 files changed, 150 insertions(+), 84 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/controllers/cross.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/controllers/cross.go	Mon Oct 14 14:58:04 2019 +0200
@@ -78,7 +78,7 @@
 		ctx, conn,
 		csi.Properties.Bottleneck, csi.Properties.Date.Time)
 
-	log.Printf("info: loading octree took %s\n", time.Since(start))
+	log.Printf("info: loading mesh took %s\n", time.Since(start))
 	if err != nil {
 		return
 	}
@@ -93,7 +93,7 @@
 		return
 	}
 
-	// The coordinate system of the octree is an UTM projection.
+	// The coordinate system of the mesh is an UTM projection.
 	// The input coordinates are in WGS84.
 	// So we need to reproject them.
 
@@ -102,7 +102,7 @@
 	var rp *models.Reprojector
 	if rp, err = models.NewReprojector(
 		ctx, conn,
-		models.WGS84, tree.EPSG,
+		models.WGS84, tree.EPSG(),
 	); err != nil {
 		return
 	}
@@ -142,14 +142,14 @@
 		}
 
 	}
-	log.Printf("info: octree traversal took %s\n", time.Since(start))
+	log.Printf("info: mesh traversal took %s\n", time.Since(start))
 
 	start = time.Now()
 
 	var joined models.GeoJSONMultiLineCoordinatesZ
 	joined, err = projectBack(
 		ctx,
-		segments, tree.EPSG,
+		segments, tree.EPSG(),
 		conn,
 	)
 
--- a/pkg/controllers/diff.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/controllers/diff.go	Mon Oct 14 14:58:04 2019 +0200
@@ -134,7 +134,7 @@
 		ctx, conn,
 		dci.Bottleneck, dci.Minuend.Time)
 
-	log.Printf("info: loading minuend octree took %s\n", time.Since(start))
+	log.Printf("info: loading minuend mesh took %s\n", time.Since(start))
 	if err != nil {
 		return
 	}
@@ -155,7 +155,7 @@
 		ctx, conn,
 		dci.Bottleneck, dci.Subtrahend.Time)
 
-	log.Printf("info: loading subtrahend octree took %s\n", time.Since(start))
+	log.Printf("info: loading subtrahend mesh took %s\n", time.Since(start))
 	if err != nil {
 		return
 	}
@@ -171,12 +171,12 @@
 	}
 
 	// We need a slow path implementation for this.
-	epsg := minuendTree.EPSG
-	if epsg != subtrahendTree.EPSG {
+	epsg := minuendTree.EPSG()
+	if epsg != subtrahendTree.EPSG() {
 		err = mw.JSONError{
 			Code: http.StatusInternalServerError,
 			Message: "Calculating differences between two different " +
-				"EPSG code octrees are not supported, yet.",
+				"EPSG code meshes are not supported, yet.",
 		}
 		return
 	}
@@ -254,7 +254,7 @@
 
 	tree.BuildWithout(tin, removed)
 
-	log.Printf("info: Building final STRTree took: %v\n", time.Since(start))
+	log.Printf("info: Building final mesh took: %v\n", time.Since(start))
 
 	start = time.Now()
 
--- a/pkg/imports/isr.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/imports/isr.go	Mon Oct 14 14:58:04 2019 +0200
@@ -171,11 +171,11 @@
 
 	// For all sounding results in bottleneck.
 	for _, sr := range bn.srs {
-		tree, err := octree.FetchOctreeDirectly(ctx, tx, sr)
+		tree, err := octree.FetchMeshDirectly(ctx, tx, sr)
 		if err != nil {
 			return err
 		}
-		hs := octree.ExtrapolateClassBreaks(heights, tree.Min.Z, tree.Max.Z)
+		hs := octree.ExtrapolateClassBreaks(heights, tree.Min().Z, tree.Max().Z)
 		hs = common.DedupFloat64s(hs)
 
 		// Delete the old iso areas.
@@ -184,7 +184,7 @@
 		}
 
 		// Calculate and store the iso areas.
-		areas := octree.TraceAreas(hs, isoCellSize, tree.Min, tree.Max, tree.Value)
+		areas := octree.TraceAreas(hs, isoCellSize, tree.Min(), tree.Max(), tree.Value)
 		for i, a := range areas {
 			if len(a) == 0 {
 				continue
--- a/pkg/octree/cache.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/octree/cache.go	Mon Oct 14 14:58:04 2019 +0200
@@ -28,7 +28,7 @@
 
 	cacheEntry struct {
 		checksum string
-		tree     *Tree
+		tree     *STRTree
 		access   time.Time
 	}
 
@@ -47,24 +47,24 @@
 )
 
 const (
-	directFetchOctreeSQL = `
-SELECT octree_index FROM waterway.sounding_results
+	directMeshSQL = `
+SELECT mesh_index FROM waterway.sounding_results
 WHERE id = $1
 `
-	fetchOctreeSQL = `
-SELECT octree_checksum, octree_index
+	fetchMeshSQL = `
+SELECT mesh_checksum, mesh_index
 FROM waterway.sounding_results
 WHERE bottleneck_id = $1 AND date_info = $2::date
-  AND octree_checksum IS NOT NULL AND octree_index IS NOT NULL
+  AND mesh_checksum IS NOT NULL AND mesh_index IS NOT NULL
 `
-	checkOctreeSQL = `
+	checkMeshSQL = `
 SELECT CASE
-  WHEN octree_checksum = $3 THEN NULL
-  ELSE octree_index
+  WHEN mesh_checksum = $3 THEN NULL
+  ELSE mesh_index
   END
 FROM waterway.sounding_results
 WHERE bottleneck_id = $1 AND date_info = $2::date
-  AND octree_checksum IS NOT NULL AND octree_index IS NOT NULL
+  AND mesh_checksum IS NOT NULL AND mesh_index IS NOT NULL
 `
 )
 
@@ -99,29 +99,33 @@
 	ctx context.Context,
 	conn *sql.Conn,
 	bottleneck string, date time.Time,
-) (*Tree, error) {
+) (*STRTree, error) {
 	return cache.get(ctx, conn, bottleneck, date)
 }
 
 // FetchOctreeDirectly loads an octree directly from the database.
-func FetchOctreeDirectly(
+func FetchMeshDirectly(
 	ctx context.Context,
 	tx *sql.Tx,
 	id int64,
-) (*Tree, error) {
+) (*STRTree, error) {
 	var data []byte
-	err := tx.QueryRowContext(ctx, directFetchOctreeSQL, id).Scan(&data)
+	err := tx.QueryRowContext(ctx, directMeshSQL, id).Scan(&data)
 	if err != nil {
 		return nil, err
 	}
-	return Deserialize(data)
+	tree := new(STRTree)
+	if err := tree.FromBytes(data); err != nil {
+		return nil, err
+	}
+	return tree, nil
 }
 
 func (c *Cache) get(
 	ctx context.Context,
 	conn *sql.Conn,
 	bottleneck string, date time.Time,
-) (*Tree, error) {
+) (*STRTree, error) {
 	c.Lock()
 	defer c.Unlock()
 
@@ -134,7 +138,7 @@
 	if entry == nil {
 		// fetch from database
 		err := conn.QueryRowContext(
-			ctx, fetchOctreeSQL, bottleneck, date).Scan(&checksum, &data)
+			ctx, fetchMeshSQL, bottleneck, date).Scan(&checksum, &data)
 		switch {
 		case err == sql.ErrNoRows:
 			return nil, nil
@@ -144,7 +148,7 @@
 	} else {
 		// check if we are not outdated.
 		err := conn.QueryRowContext(
-			ctx, checkOctreeSQL, bottleneck, date, entry.checksum).Scan(&data)
+			ctx, checkMeshSQL, bottleneck, date, entry.checksum).Scan(&data)
 		switch {
 		case err == sql.ErrNoRows:
 			return nil, nil
@@ -157,8 +161,9 @@
 		}
 	}
 
-	tree, err := Deserialize(data)
-	if err != nil {
+	tree := new(STRTree)
+
+	if err := tree.FromBytes(data); err != nil {
 		return nil, err
 	}
 
--- a/pkg/octree/loader.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/octree/loader.go	Mon Oct 14 14:58:04 2019 +0200
@@ -29,6 +29,7 @@
 		return err
 	}
 	index := make([]int32, numIndex)
+	s.index = index
 
 	var last int32
 	for i := range index {
--- a/pkg/octree/strtree.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/octree/strtree.go	Mon Oct 14 14:58:04 2019 +0200
@@ -32,12 +32,72 @@
 	bboxes  []Box2D
 }
 
+// Vertical does a vertical cross cut from (x1, y1) to (x2, y2).
+func (s *STRTree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {
+	box := Box2D{
+		X1: math.Min(x1, x2),
+		Y1: math.Min(y1, y2),
+		X2: math.Max(x1, x2),
+		Y2: math.Max(y1, y2),
+	}
+
+	// out of bounding box
+	if box.X2 < s.tin.Min.X || s.tin.Max.X < box.X1 ||
+		box.Y2 < s.tin.Min.Y || s.tin.Max.Y < box.Y1 {
+		return
+	}
+
+	line := NewPlane2D(x1, y1, x2, y2)
+
+	vertices := s.tin.Vertices
+
+	stack := []int32{s.index[0]}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		if top > 0 { // node
+			if b := s.bbox(top); b.Intersects(box) && b.IntersectsPlane(line) {
+				n := s.index[top+1]
+				stack = append(stack, s.index[top+2:top+2+n]...)
+			}
+		} else { // leaf
+			top = -top - 1
+			if b := s.bbox(top); !b.Intersects(box) || !b.IntersectsPlane(line) {
+				continue
+			}
+			n := s.index[top+1]
+			for _, idx := range s.index[top+2 : top+2+n] {
+				ti := s.tin.Triangles[idx]
+				t := Triangle{
+					vertices[ti[0]],
+					vertices[ti[1]],
+					vertices[ti[2]],
+				}
+				v0 := line.Eval(t[0].X, t[0].Y)
+				v1 := line.Eval(t[1].X, t[1].Y)
+				v2 := line.Eval(t[2].X, t[2].Y)
+
+				if onPlane(v0) || onPlane(v1) || onPlane(v2) ||
+					sides(sides(sides(0, v0), v1), v2) == 3 {
+					fn(&t)
+				}
+			}
+		}
+	}
+}
+
 func (s *STRTree) Min() Vertex  { return s.tin.Min }
 func (s *STRTree) Max() Vertex  { return s.tin.Max }
 func (s *STRTree) EPSG() uint32 { return s.tin.EPSG }
 
 func (s *STRTree) Value(x, y float64) (float64, bool) {
 
+	if len(s.index) == 0 {
+		return 0, false
+	}
+
 	stack := []int32{s.index[0]}
 
 	vertices := s.tin.Vertices
@@ -56,8 +116,8 @@
 			if !s.bbox(top).Contains(x, y) {
 				continue
 			}
-			for i, n := int32(0), s.index[top+1]; i < n; i++ {
-				idx := s.index[top+2+i]
+			n := s.index[top+1]
+			for _, idx := range s.index[top+2 : top+2+n] {
 				ti := s.tin.Triangles[idx]
 				t := Triangle{
 					vertices[ti[0]],
@@ -412,3 +472,51 @@
 	}
 	return int32(-(pos + 1))
 }
+
+func (s *STRTree) Diff(other *STRTree) PointMap {
+
+	firstVs, secondVs := s.tin.Vertices, other.tin.Vertices
+
+	result := make(PointMap, len(firstVs)+len(secondVs))
+
+	sliceWork(
+		firstVs,
+		result,
+		func(slice []Vertex, turn func([]Vertex) []Vertex) {
+			p := turn(nil)
+			for i := range slice {
+				v := &slice[i]
+				if z, found := other.Value(v.X, v.Y); found {
+					p = append(p, Vertex{v.X, v.Y, v.Z - z})
+					if len(p) == cap(p) {
+						p = turn(p)
+					}
+				}
+			}
+			if len(p) > 0 {
+				turn(p)
+			}
+		})
+
+	sliceWork(
+		secondVs,
+		result,
+		func(
+			slice []Vertex, turn func([]Vertex) []Vertex) {
+			p := turn(nil)
+			for i := range slice {
+				v := &slice[i]
+				if z, found := s.Value(v.X, v.Y); found {
+					p = append(p, Vertex{v.X, v.Y, z - v.Z})
+					if len(p) == cap(p) {
+						p = turn(p)
+					}
+				}
+			}
+			if len(p) > 0 {
+				turn(p)
+			}
+		})
+
+	return result
+}
--- a/pkg/octree/tree.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/octree/tree.go	Mon Oct 14 14:58:04 2019 +0200
@@ -312,51 +312,3 @@
 		}
 	}
 }
-
-func (ot *Tree) Diff(other *Tree) PointMap {
-
-	firstVs, secondVs := ot.Vertices(), other.Vertices()
-
-	result := make(PointMap, len(firstVs)+len(secondVs))
-
-	sliceWork(
-		firstVs,
-		result,
-		func(slice []Vertex, turn func([]Vertex) []Vertex) {
-			p := turn(nil)
-			for i := range slice {
-				v := &slice[i]
-				if z, found := other.Value(v.X, v.Y); found {
-					p = append(p, Vertex{v.X, v.Y, v.Z - z})
-					if len(p) == cap(p) {
-						p = turn(p)
-					}
-				}
-			}
-			if len(p) > 0 {
-				turn(p)
-			}
-		})
-
-	sliceWork(
-		secondVs,
-		result,
-		func(
-			slice []Vertex, turn func([]Vertex) []Vertex) {
-			p := turn(nil)
-			for i := range slice {
-				v := &slice[i]
-				if z, found := ot.Value(v.X, v.Y); found {
-					p = append(p, Vertex{v.X, v.Y, z - v.Z})
-					if len(p) == cap(p) {
-						p = turn(p)
-					}
-				}
-			}
-			if len(p) > 0 {
-				turn(p)
-			}
-		})
-
-	return result
-}