Mercurial > gemma
changeset 4663:aa3a1e7e9527
Merged stree-experiment back into default.
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Mon, 14 Oct 2019 15:58:16 +0200 |
| parents | a89e4db7980b (current diff) a2f8b3ad237a (diff) |
| children | 7d2463c7b4ad |
| files | |
| diffstat | 16 files changed, 708 insertions(+), 210 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/oct2str/main.go Mon Oct 14 15:58:16 2019 +0200 @@ -0,0 +1,167 @@ +// 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) 2019 by via donau +// – Österreichische Wasserstraßen-Gesellschaft mbH +// Software engineering by Intevation GmbH +// +// Author(s): +// * Sascha L. Teichmann <sascha.teichmann@intevation.de> + +package main + +import ( + "context" + "crypto/sha1" + "database/sql" + "encoding/hex" + "flag" + "log" + + "gemma.intevation.de/gemma/pkg/octree" + "github.com/jackc/pgx" + "github.com/jackc/pgx/stdlib" +) + +const ( + fetchOneOctreeSQL = ` +SELECT + id, + octree_index +FROM waterway.sounding_results +WHERE mesh_index IS NULL +LIMIT 1` + + storeSTRTreeSQL = ` +UPDATE waterway.sounding_results +SET mesh_index = $1, mesh_checksum = $2 +WHERE id = $3` +) + +func process(ctx context.Context, conn *sql.Conn, count int) error { + + fetch, err := conn.PrepareContext(ctx, fetchOneOctreeSQL) + if err != nil { + return err + } + defer fetch.Close() + + store, err := conn.PrepareContext(ctx, storeSTRTreeSQL) + if err != nil { + return err + } + defer store.Close() + + var next func() bool + if count < 0 { + next = func() bool { return true } + } else { + var c int + next = func() bool { + if c < count { + c++ + return true + } + return false + } + } + +loop: + for next() { + switch err := func() error { + tx, err := conn.BeginTx(ctx, nil) + if err != nil { + return err + } + defer tx.Rollback() + + var id int64 + var data []byte + + if err = tx.Stmt(fetch).QueryRowContext(ctx).Scan(&id, &data); err != nil { + return err + } + + otree, err := octree.Deserialize(data) + if err != nil { + return err + } + + unused := otree.FindUnused() + + log.Printf("unused: %d\n", len(unused)) + + str := octree.STRTree{Entries: 16} + str.BuildWithout(otree.Tin(), unused) + + out, err := str.Bytes() + if err != nil { + return err + } + h := sha1.New() + h.Write(out) + checksum := hex.EncodeToString(h.Sum(nil)) + log.Printf("hash: %s\n", checksum) + + if _, err := tx.Stmt(store).ExecContext(ctx, out, checksum, id); err != nil { + return err + } + + return tx.Commit() + }(); { + case err == sql.ErrNoRows: + break loop + case err != nil: + return err + } + } + + return nil +} + +func connect(cc pgx.ConnConfig, count int) error { + db := stdlib.OpenDB(cc) + defer db.Close() + + ctx := context.Background() + + conn, err := db.Conn(ctx) + if err != nil { + return err + } + defer conn.Close() + + return process(ctx, conn, count) +} + +func main() { + var ( + db = flag.String("database", "gemma", "database name") + user = flag.String("user", "sophie", "database user") + host = flag.String("host", "localhost", "database host") + password = flag.String("password", "so2Phie4", "database password") + port = flag.Uint("port", 5432, "database port") + ssl = flag.String("ssl", "prefer", "SSL mode") + count = flag.Int("count", -1, "how many sounding results to convert") + ) + + flag.Parse() + cc, err := pgx.ParseConnectionString("sslmode=" + *ssl) + if err != nil { + log.Fatalf("error: %v\n", err) + } + + // Do the rest manually to allow whitespace in user/password. + cc.Host = *host + cc.Port = uint16(*port) + cc.User = *user + cc.Password = *password + cc.Database = *db + + if err := connect(cc, *count); err != nil { + log.Fatalf("error: %v\n", err) + } +}
--- a/pkg/controllers/cross.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/controllers/cross.go Mon Oct 14 15:58:16 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 14:44:51 2019 +0200 +++ b/pkg/controllers/diff.go Mon Oct 14 15:58:16 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 14:44:51 2019 +0200 +++ b/pkg/imports/isr.go Mon Oct 14 15:58:16 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/imports/sr.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/imports/sr.go Mon Oct 14 15:58:16 2019 +0200 @@ -153,9 +153,9 @@ ST_AsBinary(ST_Buffer(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer), 0.1)), ST_Area(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))` - insertOctreeSQL = ` + insertMeshSQL = ` UPDATE waterway.sounding_results SET - octree_checksum = $2, octree_index = $3 + mesh_checksum = $2, mesh_index = $3 WHERE id = $1` repairBoundarySQL = ` @@ -497,12 +497,9 @@ start = time.Now() - var tree *octree.Tree - { - builder := octree.NewBuilder(tri.Tin()) - builder.Build(removed) - tree = builder.Tree() - } + tin := tri.Tin() + var virtual octree.STRTree + virtual.BuildWithout(tin, removed) feedback.Info("Building took %v", time.Since(start)) @@ -517,9 +514,13 @@ generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) - tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) { - generated = append(generated, vertices...) - }) + octree.GenerateRandomVertices( + numPoints, + tin.Min, tin.Max, + virtual.Value, + func(vertices []octree.Vertex) { + generated = append(generated, vertices...) + }) feedback.Info("Generating %d points took %v.", len(generated), time.Since(start)) @@ -531,7 +532,7 @@ return } dupes[key] = struct{}{} - if z, ok := tree.Value(x, y); ok { + if z, ok := virtual.Value(x, y); ok { generated = append(generated, octree.Vertex{X: x, Y: y, Z: z}) } }) @@ -547,7 +548,6 @@ feedback.Info("Second triangulation took %v.", time.Since(start)) feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) feedback.Info("Clipping triangles from new mesh.") - } start = time.Now() @@ -556,25 +556,21 @@ var str octree.STRTree str.Build(tin) - feedback.Info("Building STR tree took %v", time.Since(start)) + feedback.Info("Building clipping index took %v", time.Since(start)) start = time.Now() clippingPolygonBuffered.Indexify() removed = str.Clip(&clippingPolygonBuffered) - feedback.Info("Clipping STR tree took %v.", time.Since(start)) + feedback.Info("Clipping took %v.", time.Since(start)) feedback.Info("Number of triangles to clip %d.", len(removed)) - feedback.Info("Build final octree index") + start = time.Now() + final := octree.STRTree{Entries: 16} + final.BuildWithout(tin, removed) - builder := octree.NewBuilder(tin) - builder.Build(removed) - octreeIndex, err := builder.Bytes() - if err != nil { - return nil, err - } - feedback.Info("Building octree took %v.", time.Since(start)) - feedback.Info("Store octree.") + feedback.Info("Building final mesh took %v.", time.Since(start)) + feedback.Info("Store mesh.") start = time.Now() @@ -612,15 +608,20 @@ return nil, err } - h := sha1.New() - h.Write(octreeIndex) - checksum := hex.EncodeToString(h.Sum(nil)) - _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex) + index, err := final.Bytes() if err != nil { return nil, err } - feedback.Info("Storing octree index took %s.", time.Since(start)) - err = generateIsos(ctx, tx, feedback, builder.Tree(), id) + + h := sha1.New() + h.Write(index) + checksum := hex.EncodeToString(h.Sum(nil)) + _, err = tx.ExecContext(ctx, insertMeshSQL, id, checksum, index) + if err != nil { + return nil, err + } + feedback.Info("Storing mesh index took %s.", time.Since(start)) + err = generateIsos(ctx, tx, feedback, &final, id) if err != nil { return nil, err } @@ -827,7 +828,7 @@ ctx context.Context, tx *sql.Tx, feedback Feedback, - tree *octree.Tree, + tree *octree.STRTree, id int64, ) error { @@ -836,16 +837,18 @@ "morphology_classbreaks", ) + minZ, maxZ := tree.Min().Z, tree.Max().Z + if err != nil { feedback.Warn("Loading class breaks failed: %v", err) feedback.Info("Using default class breaks") heights = nil - h := contourStepWidth * math.Ceil(tree.Min.Z/contourStepWidth) - for ; h <= tree.Max.Z; h += contourStepWidth { + h := contourStepWidth * math.Ceil(minZ/contourStepWidth) + for ; h <= maxZ; h += contourStepWidth { heights = append(heights, h) } } else { - heights = octree.ExtrapolateClassBreaks(heights, tree.Min.Z, tree.Max.Z) + heights = octree.ExtrapolateClassBreaks(heights, minZ, maxZ) // We set steps for InBetweenClassBreaks to 1, so it // becomes a null operation. The extra class breaks // were considered unexpected and confusing by the @@ -870,7 +873,7 @@ ctx context.Context, tx *sql.Tx, feedback Feedback, - tree *octree.Tree, + tree *octree.STRTree, heights []float64, id int64, ) error { @@ -881,11 +884,11 @@ time.Since(total)) }() - areas := octree.TraceAreas(heights, isoCellSize, tree.Min, tree.Max, tree.Value) + areas := octree.TraceAreas(heights, isoCellSize, tree.Min(), tree.Max(), tree.Value) return storeAreas( ctx, tx, feedback, - areas, tree.EPSG, heights, id) + areas, tree.EPSG(), heights, id) } func storeAreas(
--- a/pkg/octree/areas.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/octree/areas.go Mon Oct 14 15:58:16 2019 +0200 @@ -26,6 +26,73 @@ "gemma.intevation.de/gemma/pkg/wkb" ) +func GenerateRandomVertices( + n int, + min, max Vertex, + eval func(float64, float64) (float64, bool), + callback func([]Vertex), +) { + var wg sync.WaitGroup + + jobs := make(chan int) + out := make(chan []Vertex) + done := make(chan struct{}) + + cpus := runtime.NumCPU() + + free := make(chan []Vertex, cpus) + + for i := 0; i < cpus; i++ { + wg.Add(1) + go func() { + defer wg.Done() + xRange := common.Random(min.X, max.X) + yRange := common.Random(min.Y, max.Y) + + for size := range jobs { + var vertices []Vertex + select { + case vertices = <-free: + default: + vertices = make([]Vertex, 0, 1000) + } + for len(vertices) < size { + x, y := xRange(), yRange() + if z, ok := eval(x, y); ok { + vertices = append(vertices, Vertex{X: x, Y: y, Z: z}) + } + } + out <- vertices + } + }() + } + + go func() { + defer close(done) + for vertices := range out { + callback(vertices) + select { + case free <- vertices[:0]: + default: + } + } + }() + + for remain := n; remain > 0; { + if remain > 1000 { + jobs <- 1000 + remain -= 1000 + } else { + jobs <- remain + remain = 0 + } + } + close(jobs) + wg.Wait() + close(out) + <-done +} + func TraceAreas( heights []float64, cellSize float64,
--- a/pkg/octree/cache.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/octree/cache.go Mon Oct 14 15:58:16 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 14:44:51 2019 +0200 +++ b/pkg/octree/loader.go Mon Oct 14 15:58:16 2019 +0200 @@ -16,59 +16,136 @@ import ( "bufio" "bytes" + "compress/gzip" "encoding/binary" "log" "github.com/golang/snappy" ) -func loadReader(r *bufio.Reader) (*Tree, error) { - tree := new(Tree) +func (s *STRTree) deserializeIndex(r *bufio.Reader) error { + var numIndex int32 + if err := binary.Read(r, binary.LittleEndian, &numIndex); err != nil { + return err + } + index := make([]int32, numIndex) + s.index = index + + var last int32 + for i := range index { + v, err := binary.ReadVarint(r) + if err != nil { + return err + } + value := int32(v) + last + index[i] = value + last = value + } - if err := binary.Read(r, binary.LittleEndian, &tree.EPSG); err != nil { - return nil, err + return nil +} + +func (s *STRTree) deserializeBBoxes(r *bufio.Reader) error { + + var numBBoxes int32 + if err := binary.Read(r, binary.LittleEndian, &numBBoxes); err != nil { + return err + } + + bboxes := make([]Box2D, numBBoxes) + s.bboxes = bboxes + + var err error + + read := func(v *float64) { + if err == nil { + err = binary.Read(r, binary.LittleEndian, v) + } } - log.Printf("info: EPSG: %d\n", tree.EPSG) + for i := range bboxes { + read(&bboxes[i].X1) + read(&bboxes[i].Y1) + read(&bboxes[i].X2) + read(&bboxes[i].Y2) + } + + return err +} - if err := tree.Min.Read(r); err != nil { - return nil, err +func (s *STRTree) deserialize(r *bufio.Reader) error { + s.tin = new(Tin) + + if err := s.tin.Deserialize(r); err != nil { + return err + } + var numEntries uint8 + if err := binary.Read(r, binary.LittleEndian, &numEntries); err != nil { + return err + } + s.Entries = int(numEntries) + + if err := s.deserializeIndex(r); err != nil { + return err } - if err := tree.Max.Read(r); err != nil { - return nil, err + return s.deserializeBBoxes(r) +} + +func (s *STRTree) FromBytes(data []byte) error { + r, err := gzip.NewReader(bytes.NewReader(data)) + if err != nil { + return err + } + return s.deserialize(bufio.NewReader(r)) +} + +func (t *Tin) Deserialize(r *bufio.Reader) error { + + if err := binary.Read(r, binary.LittleEndian, &t.EPSG); err != nil { + return err + } + + log.Printf("info: EPSG: %d\n", t.EPSG) + + if err := t.Min.Read(r); err != nil { + return err + } + + if err := t.Max.Read(r); err != nil { + return err } log.Printf("info: BBOX: [[%f, %f, %f], [%f, %f, %f]]\n", - tree.Min.X, tree.Min.Y, tree.Min.Z, - tree.Max.X, tree.Max.Y, tree.Max.Z) + t.Min.X, t.Min.Y, t.Min.Z, + t.Max.X, t.Max.Y, t.Max.Z) var numVertices uint32 if err := binary.Read(r, binary.LittleEndian, &numVertices); err != nil { - return nil, err + return err } log.Printf("info: vertices: %d\n", numVertices) vertices := make([]Vertex, numVertices) - tree.vertices = vertices + t.Vertices = vertices for i := range vertices { if err := vertices[i].Read(r); err != nil { - return nil, err + return err } } var numTriangles uint32 if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil { - return nil, err + return err } log.Printf("info: triangles: %d\n", numTriangles) indices := make([]int32, 3*numTriangles) triangles := make([][]int32, numTriangles) - tree.triangles = triangles + t.Triangles = triangles var last int32 @@ -79,7 +156,7 @@ for j := range tri { v, err := binary.ReadVarint(r) if err != nil { - return nil, err + return err } value := int32(v) + last tri[j] = value @@ -87,6 +164,24 @@ } } + return nil +} + +func loadReader(r *bufio.Reader) (*Tree, error) { + tree := new(Tree) + + var tin Tin + + if err := tin.Deserialize(r); err != nil { + return nil, err + } + + tree.EPSG = tin.EPSG + tree.vertices = tin.Vertices + tree.triangles = tin.Triangles + tree.Min = tin.Min + tree.Max = tin.Max + var numNodes uint32 if err := binary.Read(r, binary.LittleEndian, &numNodes); err != nil { return nil, err @@ -97,7 +192,7 @@ tree.index = make([]int32, numNodes) entries := tree.index[1:] - last = 0 + var last int32 for i := range entries { v, err := binary.ReadVarint(r) if err != nil {
--- a/pkg/octree/strtree.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/octree/strtree.go Mon Oct 14 15:58:16 2019 +0200 @@ -14,6 +14,11 @@ package octree import ( + "bytes" + "compress/gzip" + "encoding/binary" + "io" + "log" "math" "sort" ) @@ -27,8 +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 @@ -47,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]], @@ -153,6 +222,94 @@ return removed } +func (s *STRTree) serializeIndex(w io.Writer) error { + + if err := binary.Write(w, binary.LittleEndian, int32(len(s.index))); err != nil { + return err + } + + var buf [binary.MaxVarintLen32]byte + + var last int32 + var written int + + for _, x := range s.index { + delta := x - last + n := binary.PutVarint(buf[:], int64(delta)) + for p := buf[:n]; len(p) > 0; p = p[n:] { + var err error + if n, err = w.Write(p); err != nil { + return err + } + written += n + } + last = x + } + + log.Printf("info: compressed index in bytes: %d %.2f (%d %.2f)\n", + written, + float64(written)/(1024*1024), + 4*len(s.index), + float64(4*len(s.index))/(1024*1024), + ) + + return nil +} + +func (s *STRTree) serializeBBoxes(w io.Writer) (rr error) { + + if err := binary.Write(w, binary.LittleEndian, int32(len(s.bboxes))); err != nil { + return err + } + + var err error + + write := func(v float64) { + if err == nil { + err = binary.Write(w, binary.LittleEndian, math.Float64bits(v)) + } + } + for _, box := range s.bboxes { + write(box.X1) + write(box.Y1) + write(box.X2) + write(box.Y2) + } + + return err +} + +func (s *STRTree) Bytes() ([]byte, error) { + + var buf bytes.Buffer + w, err := gzip.NewWriterLevel(&buf, gzip.BestSpeed) + if err != nil { + return nil, err + } + + if err := s.tin.serialize(w); err != nil { + return nil, err + } + + if err := binary.Write(w, binary.LittleEndian, uint8(s.Entries)); err != nil { + return nil, err + } + + if err := s.serializeIndex(w); err != nil { + return nil, err + } + + if err := s.serializeBBoxes(w); err != nil { + return nil, err + } + + if err := w.Close(); err != nil { + return nil, err + } + + return buf.Bytes(), nil +} + func (s *STRTree) allTriangles(pos int32, tris map[int32]struct{}) { stack := []int32{pos} for len(stack) > 0 { @@ -315,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/tin.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/octree/tin.go Mon Oct 14 15:58:16 2019 +0200 @@ -224,10 +224,21 @@ return err } + var err error + vwrite := func(v float64) { + if err == nil { + err = binary.Write(w, binary.LittleEndian, math.Float64bits(v)) + } + } + for _, v := range t.Vertices { - if err := v.Write(w); err != nil { - return err - } + vwrite(v.X) + vwrite(v.Y) + vwrite(v.Z) + } + + if err != nil { + return err } log.Printf("info: vertices %d (%d)\n", len(t.Vertices), len(t.Vertices)*3*8)
--- a/pkg/octree/tree.go Mon Oct 14 14:44:51 2019 +0200 +++ b/pkg/octree/tree.go Mon Oct 14 15:58:16 2019 +0200 @@ -15,10 +15,6 @@ import ( "math" - "runtime" - "sync" - - "gemma.intevation.de/gemma/pkg/common" ) // Tree is an Octree holding triangles. @@ -52,6 +48,50 @@ {0.5, 0.5, 1.0, 1.0}, } +func (ot *Tree) Tin() *Tin { + return &Tin{ + EPSG: ot.EPSG, + Vertices: ot.vertices, + Triangles: ot.triangles, + Min: ot.Min, + Max: ot.Max, + } +} + +func (ot *Tree) FindUnused() map[int32]struct{} { + + used := make(map[int32]struct{}, len(ot.triangles)) + for i := int32(0); i < int32(len(ot.triangles)); i++ { + used[i] = struct{}{} + } + + stack := []int32{1} + + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + + if top > 0 { // node + if index := ot.index[top:]; len(index) > 7 { + for _, idx := range index[:8] { + if idx != 0 { + stack = append(stack, idx) + } + } + } + } else { // leaf + pos := -top - 1 + n := ot.index[pos] + indices := ot.index[pos+1 : pos+1+n] + for _, idx := range indices { + delete(used, idx) + } + } + } + + return used +} + func (ot *Tree) Value(x, y float64) (float64, bool) { // out of bounding box @@ -272,113 +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 -} - -func (ot *Tree) GenerateRandomVertices(n int, callback func([]Vertex)) { - var wg sync.WaitGroup - - jobs := make(chan int) - out := make(chan []Vertex) - done := make(chan struct{}) - - cpus := runtime.NumCPU() - - free := make(chan []Vertex, cpus) - - for i := 0; i < cpus; i++ { - wg.Add(1) - go func() { - defer wg.Done() - xRange := common.Random(ot.Min.X, ot.Max.X) - yRange := common.Random(ot.Min.Y, ot.Max.Y) - - for size := range jobs { - var vertices []Vertex - select { - case vertices = <-free: - default: - vertices = make([]Vertex, 0, 1000) - } - for len(vertices) < size { - x, y := xRange(), yRange() - if z, ok := ot.Value(x, y); ok { - vertices = append(vertices, Vertex{X: x, Y: y, Z: z}) - } - } - out <- vertices - } - }() - } - - go func() { - defer close(done) - for vertices := range out { - callback(vertices) - select { - case free <- vertices[:0]: - default: - } - } - }() - - for remain := n; remain > 0; { - if remain > 1000 { - jobs <- 1000 - remain -= 1000 - } else { - jobs <- remain - remain = 0 - } - } - close(jobs) - wg.Wait() - close(out) - <-done -}
--- a/schema/gemma.sql Mon Oct 14 14:44:51 2019 +0200 +++ b/schema/gemma.sql Mon Oct 14 15:58:16 2019 +0200 @@ -712,6 +712,8 @@ depth_reference varchar NOT NULL, -- REFERENCES depth_references, octree_checksum varchar, octree_index bytea, + mesh_checksum varchar, + mesh_index bytea, staging_done boolean NOT NULL DEFAULT false ) CREATE CONSTRAINT TRIGGER a_sounding_results_reference_bottleneck
--- a/schema/update-db.sh Mon Oct 14 14:44:51 2019 +0200 +++ b/schema/update-db.sh Mon Oct 14 15:58:16 2019 +0200 @@ -104,7 +104,7 @@ echo "Running updates for $new_ver ..." file_args="" - for f in "$d"/* ; do + for f in "$d"/*.sql ; do file_args+=" -f $f" done
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/updates/1302/01.mesh-columns.sql Mon Oct 14 15:58:16 2019 +0200 @@ -0,0 +1,2 @@ +alter table waterway.sounding_results add column mesh_index bytea; +alter table waterway.sounding_results add column mesh_checksum varchar;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/updates/1302/after-update.sql_manual Mon Oct 14 15:58:16 2019 +0200 @@ -0,0 +1,11 @@ +-- +-- Before executing this script you have to build the cmd/oct2str tool. +-- $ cd cmd/oct2str && go build +-- Use +-- $ ./oct2str --help +-- to find out how to use it: +-- Passing the right db credentials should be enough. +-- + +ALTER TABLE waterway.sounding_results DROP COLUMN octree_index; +ALTER TABLE waterway.sounding_results DROP COLUMN octree_checksum;