# HG changeset patch # User Sascha L. Teichmann # Date 1571057884 -7200 # Node ID 4bbfe3dd2ab5d8a94ceb90e119d0f4ccc0222cad # Parent a2f16987911bed13943870e0db0cec54ccfc9f01 Completed usage of STRTrees. diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/controllers/cross.go --- 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, ) diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/controllers/diff.go --- 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() diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/imports/isr.go --- 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 diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/octree/cache.go --- 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 } diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/octree/loader.go --- 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 { diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/octree/strtree.go --- 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 +} diff -r a2f16987911b -r 4bbfe3dd2ab5 pkg/octree/tree.go --- 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 -}