Mercurial > gemma
changeset 2481:3cf5d27a6c8b octree-diff
Merged defualt into octree-diff branch.
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Fri, 01 Mar 2019 11:06:27 +0100 |
| parents | 242104c338ff (diff) 9de710bdb535 (current diff) |
| children | 620038ade708 |
| files | client/src/components/importqueue/Importqueue.vue client/src/components/importqueue/Importqueuedetail.vue client/src/components/staging/Staging.vue client/src/components/staging/StagingDetail.vue client/src/components/ui/box/Header.vue |
| diffstat | 13 files changed, 903 insertions(+), 97 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octreediff/db.go Fri Mar 01 11:06:27 2019 +0100 @@ -0,0 +1,52 @@ +// 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 main + +import ( + "database/sql" + "flag" + + "github.com/jackc/pgx" + "github.com/jackc/pgx/stdlib" +) + +var ( + dbhost = flag.String("dbhost", "localhost", "database host") + dbport = flag.Uint("dbport", 5432, "database port") + dbname = flag.String("dbname", "gemma", "database user") + dbuser = flag.String("dbuser", "scott", "database user") + dbpassword = flag.String("dbpw", "tiger", "database password") + dbssl = flag.String("dbssl", "prefer", "database SSL mode") +) + +func run(fn func(*sql.DB) error) error { + + // To ease SSL config ride a bit on parsing. + cc, err := pgx.ParseConnectionString("sslmode=" + *dbssl) + if err != nil { + return err + } + + // Do the rest manually to allow whitespace in user/password. + cc.Host = *dbhost + cc.Port = uint16(*dbport) + cc.User = *dbuser + cc.Password = *dbpassword + cc.Database = *dbname + + db := stdlib.OpenDB(cc) + defer db.Close() + + return fn(db) +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octreediff/main.go Fri Mar 01 11:06:27 2019 +0100 @@ -0,0 +1,471 @@ +// 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 main + +import ( + "bytes" + "context" + "database/sql" + "encoding/binary" + "errors" + "flag" + "fmt" + "log" + "math" + "runtime" + "sync" + "time" + + "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/octree" + "gemma.intevation.de/gemma/pkg/wkb" +) + +var ( + bottleneck = flag.String("bottleneck", "", "name of the bottleneck") + first = flag.String("first", "", "date of the first sounding result") + second = flag.String("second", "", "date of the second sounding result") +) + +const contourTolerance = 0.1 + +const ( + loadOctreeSQL = ` +SELECT sr.octree_index +FROM waterway.sounding_results sr JOIN waterway.bottlenecks bn + ON sr.bottleneck_id = bn.id +WHERE bn.bottleneck_id = $1 AND sr.date_info = $2::date + AND sr.octree_index IS NOT NULL` + + triangulateSQL = ` +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 +), +inter AS ( + SELECT + 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 +), +triangles AS ( + SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM ( + SELECT (ST_Dump( + ST_DelaunayTriangles(ST_GeomFromWKB($5, $2::int), 0, 2))).geom + ) t +) +SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, inter +WHERE ST_Covers(inter.area, triangles.poly) +` + + //INSERT INTO redis_diff_countour_lines ( + //INSERT INTO diff_contour_lines ( + //INSERT INTO diff_contour_lines_extern ( + insertContourSQL = ` +INSERT INTO diff_contour_lines ( + height, + lines +) +SELECT + $1, + ST_Transform( + ST_Multi( + ST_CollectionExtract( + ST_SimplifyPreserveTopology( + ST_Multi(ST_Collectionextract( + ST_MakeValid(ST_GeomFromWKB($2, $3::integer)), 2)), + $4 + ), + 2 + ) + ), + 4326 + )` +) + +func check(err error) { + if err != nil { + log.Fatalf("error: %v\n", err) + } +} + +type point struct { + x float64 + y float64 +} + +type pointMap map[point]float64 + +func (pm pointMap) asWKB() []byte { + size := 1 + 4 + 4 + len(pm)*(1+4+3*8) + + buf := bytes.NewBuffer(make([]byte, 0, size)) + + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.MultiPointZ) + binary.Write(buf, binary.LittleEndian, uint32(len(pm))) + + perPoint := bytes.NewBuffer(make([]byte, 0, 1+4)) + binary.Write(perPoint, binary.LittleEndian, wkb.NDR) + binary.Write(perPoint, binary.LittleEndian, wkb.PointZ) + hdr := perPoint.Bytes() + + for p, z := range pm { + buf.Write(hdr) + binary.Write(buf, binary.LittleEndian, math.Float64bits(p.x)) + binary.Write(buf, binary.LittleEndian, math.Float64bits(p.y)) + binary.Write(buf, binary.LittleEndian, math.Float64bits(z)) + } + + return buf.Bytes() +} + +func sliceWork( + vs []octree.Vertex, + dst pointMap, + fn func([]octree.Vertex, func([]octree.Vertex) []octree.Vertex), +) { + n := runtime.NumCPU() + + wg := new(sync.WaitGroup) + + slices := make(chan []octree.Vertex) + out := make(chan []octree.Vertex) + + pool := make(chan []octree.Vertex, n) + + const pageSize = 2048 + + turn := func(p []octree.Vertex) []octree.Vertex { + if p != nil { + out <- p + } + select { + case p = <-pool: + default: + p = make([]octree.Vertex, 0, pageSize) + } + return p + } + + for i := 0; i < n; i++ { + wg.Add(1) + go func() { + defer wg.Done() + for slice := range slices { + fn(slice, turn) + } + }() + } + done := make(chan struct{}) + go func() { + defer close(done) + for s := range out { + for i := range s { + v := &s[i] + key := point{x: v.X, y: v.Y} + if z, found := dst[key]; found { + dst[key] = (z + v.Z) * 0.5 + } else { + dst[key] = v.Z + } + } + select { + case pool <- s[:0:pageSize]: + default: + } + } + }() + + size := len(vs)/n + 1 + for len(vs) > 0 { + var l int + if len(vs) < size { + l = len(vs) + } else { + l = size + } + slices <- vs[:l] + vs = vs[l:] + } + close(slices) + wg.Wait() + close(out) + <-done +} + +func process(bottleneck string, firstDate, secondDate time.Time) error { + start := time.Now() + defer func() { log.Printf("processing took %v\n", time.Since(start)) }() + + ctx := context.Background() + + return run(func(db *sql.DB) error { + + conn, err := db.Conn(ctx) + if err != nil { + return err + } + defer conn.Close() + + tx, err := conn.BeginTx(ctx, nil) + if err != nil { + return err + } + + type load struct { + date time.Time + data []byte + err *error + dst **octree.Tree + } + + out := make(chan *load) + wg := new(sync.WaitGroup) + + n := runtime.NumCPU() + if n > 2 { + n = 2 + } + + for i := 0; i < n; i++ { + wg.Add(1) + go func() { + defer wg.Done() + for l := range out { + if *l.err == nil { + *l.dst, *l.err = octree.Deserialize(l.data) + } + } + }() + } + + var firstErr, secondErr error + var first, second *octree.Tree + + for _, l := range []*load{ + {date: firstDate, dst: &first, err: &firstErr}, + {date: secondDate, dst: &second, err: &secondErr}, + } { + var data []byte + if err := tx.QueryRowContext( + ctx, + loadOctreeSQL, + bottleneck, + l.date, + ).Scan(&data); err != nil { + *l.err = err + } else { + l.data = data + } + out <- l + } + close(out) + + wg.Wait() + + if firstErr != nil || secondErr != nil { + if firstErr != nil && secondErr != nil { + return fmt.Errorf("%v, %v", firstErr, secondErr) + } + if firstErr != nil { + return firstErr + } + return secondErr + } + + if first.EPSG != second.EPSG { + return errors.New("EPSG codes mismatch. Needs transformation slow pass.") + } + + now := time.Now() + log.Printf("loading took %v\n", now.Sub(start)) + last := now + + firstVs, secondVs := first.Vertices(), second.Vertices() + + result := make(pointMap, len(firstVs)+len(secondVs)) + + sliceWork( + firstVs, + result, + func( + slice []octree.Vertex, + turn func([]octree.Vertex) []octree.Vertex, + ) { + p := turn(nil) + for i := range slice { + v := &slice[i] + if z, found := second.Value(v.X, v.Y); found { + p = append(p, octree.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 []octree.Vertex, + turn func([]octree.Vertex) []octree.Vertex, + ) { + p := turn(nil) + for i := range slice { + v := &slice[i] + if z, found := first.Value(v.X, v.Y); found { + p = append(p, octree.Vertex{v.X, v.Y, z - v.Z}) + if len(p) == cap(p) { + p = turn(p) + } + } + } + if len(p) > 0 { + turn(p) + } + }) + + now = time.Now() + log.Printf("setting in took %v\n", now.Sub(last)) + last = now + log.Printf("num points: %d\n", len(result)) + + data := result.asWKB() + + now = time.Now() + log.Printf("turing into WKB took %v\n", now.Sub(last)) + last = now + + log.Printf("WKB size %.3fMB\n", float64(len(data))/(1024*1024)) + + var tin octree.Tin + + if err := tx.QueryRowContext( + ctx, + triangulateSQL, + bottleneck, + first.EPSG, + firstDate, + secondDate, + data, + ).Scan(&tin); err != nil { + return err + } + + now = time.Now() + log.Printf("triangulation took %v\n", now.Sub(last)) + last = now + + builder := octree.NewBuilder(&tin) + builder.Build() + + now = time.Now() + log.Printf("building octree took %v\n", now.Sub(last)) + last = now + + tree := builder.Tree() + + log.Printf("min/max: %f %f\n", tree.Min.Z, tree.Max.Z) + + var heights []float64 + + switch { + case tree.Min.Z >= 0: // All values positive. + for v := 0.0; v <= tree.Max.Z; v += 0.1 { + if v >= tree.Min.Z { + heights = append(heights, v) + } + } + case tree.Max.Z <= 0: // All values negative. + for v := 0.0; v >= tree.Min.Z; v -= 0.1 { + if v <= tree.Max.Z { + heights = append(heights, v) + } + } + default: // Positive and negative. + for v := 0.1; v <= tree.Max.Z; v += 0.1 { + heights = append(heights, v) + } + for i, j := 0, len(heights)-1; i < j; i, j = i+1, j-1 { + heights[i], heights[j] = heights[j], heights[i] + } + for v := 0.0; v >= tree.Min.Z; v -= 0.1 { + heights = append(heights, v) + } + } + + var dataSize int + + stmt, err := tx.PrepareContext(ctx, insertContourSQL) + if err != nil { + return err + } + defer stmt.Close() + + octree.DoContours(tree, heights, func(res *octree.ContourResult) { + if err == nil && len(res.Lines) > 0 { + log.Printf("%f: lines: %d\n", res.Height, len(res.Lines)) + wkb := res.Lines.AsWKB2D() + dataSize += len(wkb) + _, err = stmt.ExecContext( + ctx, + res.Height, + wkb, + first.EPSG, + contourTolerance, + ) + } + }) + + if err != nil { + return err + } + + now = time.Now() + log.Printf("Number of iso lines: %d\n", len(heights)) + log.Printf("Total WKB size: %.2fMB\n", float64(dataSize)/(1024*1024)) + log.Printf("generating iso lines took %v\n", now.Sub(last)) + last = now + + return tx.Commit() + }) +} + +func main() { + + flag.Parse() + + firstDate, err := time.Parse(common.DateFormat, *first) + check(err) + secondDate, err := time.Parse(common.DateFormat, *second) + check(err) + + if *bottleneck == "" { + log.Fatalln("Missing bottleneck name") + } + + check(process(*bottleneck, firstDate, secondDate)) +}
--- a/pkg/imports/wkb.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/imports/wkb.go Fri Mar 01 11:06:27 2019 +0100 @@ -20,6 +20,8 @@ "math" shp "github.com/jonas-p/go-shp" + + "gemma.intevation.de/gemma/pkg/wkb" ) type ( @@ -28,22 +30,14 @@ polygonSlice [][][]float64 ) -const ( - wkbNDR byte = 1 - - wkbPoint uint32 = 1 - wkbLineString uint32 = 2 - wkbPolygon uint32 = 3 -) - func (l lineSlice) asWKB() []byte { size := 1 + 4 + 4 + len(l)*(2*8) buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbLineString) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.LineString) binary.Write(buf, binary.LittleEndian, uint32(len(l))) for _, c := range l { @@ -67,8 +61,8 @@ buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbPoint) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.Point) var lat, lon float64 if len(p) > 0 { @@ -95,8 +89,8 @@ buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbPolygon) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.Polygon) binary.Write(buf, binary.LittleEndian, uint32(len(p))) for _, ring := range p {
--- a/pkg/models/cross.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/models/cross.go Fri Mar 01 11:06:27 2019 +0100 @@ -21,6 +21,8 @@ "fmt" "math" "time" + + "gemma.intevation.de/gemma/pkg/wkb" ) type ( @@ -145,22 +147,14 @@ return nil } -const ( - wkbXDR byte = 0 - wkbNDR byte = 1 - wkbLineString uint32 = 2 - wkbLineStringZ uint32 = 1000 + 2 - wkbMultiLineStringZ uint32 = 1000 + 5 -) - func (lc GeoJSONLineCoordinates) AsWKB() []byte { size := 1 + 4 + 4 + len(lc)*(2*8) buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbLineString) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.LineString) binary.Write(buf, binary.LittleEndian, uint32(len(lc))) for i := range lc { @@ -197,9 +191,9 @@ switch { case err != nil: return err - case endian == wkbNDR: + case endian == wkb.NDR: order = binary.LittleEndian - case endian == wkbXDR: + case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) @@ -211,7 +205,7 @@ switch { case err != nil: return err - case geomType != wkbLineStringZ: + case geomType != wkb.LineStringZ: return fmt.Errorf("unknown geometry type %x", geomType) } @@ -273,9 +267,9 @@ switch { case err != nil: return err - case endian == wkbNDR: + case endian == wkb.NDR: order = binary.LittleEndian - case endian == wkbXDR: + case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) @@ -287,7 +281,7 @@ switch { case err != nil: return err - case geomType != wkbMultiLineStringZ: + case geomType != wkb.MultiLineStringZ: return fmt.Errorf("unknown geometry type %d", geomType) } @@ -304,9 +298,9 @@ switch { case err != nil: return err - case endian == wkbNDR: + case endian == wkb.NDR: order = binary.LittleEndian - case endian == wkbXDR: + case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) @@ -316,7 +310,7 @@ switch { case err != nil: return err - case geomType != wkbLineStringZ: + case geomType != wkb.LineStringZ: return fmt.Errorf("unknown geometry type %d", geomType) }
--- a/pkg/octree/builder.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/builder.go Fri Mar 01 11:06:27 2019 +0100 @@ -18,6 +18,9 @@ "encoding/binary" "io" "log" + "runtime" + "sync" + "sync/atomic" "github.com/golang/snappy" ) @@ -28,8 +31,12 @@ nodes int leaves int index []int32 + + mu sync.Mutex } +type buildStep func(chan buildStep) + var cubes = [8][2]Vertex{ makeCube(0), makeCube(1), @@ -71,35 +78,154 @@ triangles[i] = int32(i) } + n := runtime.NumCPU() + + steps := make(chan buildStep) + + var wg sync.WaitGroup + for i := 0; i < n; i++ { + wg.Add(1) + go func() { + defer wg.Done() + for step := range steps { + step(steps) + } + }() + } + tb.index = append(tb.index, 0) - tb.buildRecursive(triangles, tb.t.Min, tb.t.Max, 0) + root := func(int32) { + close(steps) + } + + steps <- tb.buildConcurrent( + triangles, + tb.t.Min, tb.t.Max, + 0, + root) + + wg.Wait() + + /* + tb.buildRecursive(triangles, tb.t.Min, tb.t.Max, 0) + */ tb.index[0] = int32(len(tb.index)) log.Printf("info: num nodes: %d\n", tb.index[0]) - log.Printf("info: nodes: %d leaves: %d index %d\n", tb.nodes, tb.leaves, tb.index[0]) } +func (tb *Builder) buildConcurrent( + triangles []int32, + min, max Vertex, + depth int, + parent func(int32), +) buildStep { + + return func(steps chan buildStep) { + + // none concurrent for small parts. + if len(triangles) <= 1024 || depth > 8 { + parent(tb.buildRecursive(triangles, min, max, depth)) + return + } + + bbox := Interpolate(min, max) + + bboxes := make([][2]Vertex, len(cubes)) + + for i := range cubes { + bboxes[i] = [2]Vertex{ + bbox(cubes[i][0]), + bbox(cubes[i][1]), + } + } + + var quandrants [8][]int32 + + for _, tri := range triangles { + triangle := tb.t.Triangles[tri] + v0 := tb.t.Vertices[triangle[0]] + v1 := tb.t.Vertices[triangle[1]] + v2 := tb.t.Vertices[triangle[2]] + + l := v0 + l.Minimize(v1) + l.Minimize(v2) + + h := v0 + h.Maximize(v1) + h.Maximize(v2) + + for i := range bboxes { + if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) { + quandrants[i] = append(quandrants[i], tri) + } + } + } + + used := new(int32) + for i := range quandrants { + if len(quandrants[i]) > 0 { + *used++ + } + } + + pos := tb.allocNode() + + for i := range quandrants { + if len(quandrants[i]) > 0 { + j := int32(i) + parent := func(v int32) { + tb.index[pos+j] = v + if atomic.AddInt32(used, -1) == 0 { + parent(pos) + } + } + step := tb.buildConcurrent( + quandrants[i], + bboxes[i][0], bboxes[i][1], + depth+1, + parent) + select { + case steps <- step: + default: + // all slots busy -> execute directly. + step(steps) + } + } + } + } +} + +func (tb *Builder) allocNode() int32 { + tb.mu.Lock() + pos := int32(len(tb.index)) + tb.index = append(tb.index, + 0, 0, 0, 0, + 0, 0, 0, 0) + tb.nodes++ + tb.mu.Unlock() + return pos +} + func (tb *Builder) buildRecursive( triangles []int32, min, max Vertex, depth int, ) int32 { if len(triangles) <= 16 || depth > 8 { + tb.mu.Lock() pos := len(tb.index) tb.index = append(tb.index, int32(len(triangles))) tb.index = append(tb.index, triangles...) //log.Printf("leaf entries: %d (%d)\n", len(triangles), depth) tb.leaves++ + tb.mu.Unlock() return int32(-(pos + 1)) } - pos := len(tb.index) - tb.index = append(tb.index, - 0, 0, 0, 0, - 0, 0, 0, 0) - bbox := Interpolate(min, max) bboxes := make([][2]Vertex, len(cubes)) @@ -134,17 +260,19 @@ } } + pos := tb.allocNode() + for i := range quandrants { if len(quandrants[i]) > 0 { child := tb.buildRecursive( quandrants[i], bboxes[i][0], bboxes[i][1], depth+1) - tb.index[pos+i] = child + tb.index[pos+int32(i)] = child } } - tb.nodes++ - return int32(pos) + + return pos } func (tb *Builder) serialize(w io.Writer) error {
--- a/pkg/octree/cache.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/cache.go Fri Mar 01 11:06:27 2019 +0100 @@ -139,7 +139,7 @@ } } - tree, err := deserialize(data) + tree, err := Deserialize(data) if err != nil { return nil, err }
--- a/pkg/octree/contours.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/contours.go Fri Mar 01 11:06:27 2019 +0100 @@ -15,6 +15,7 @@ package octree import ( + "log" "runtime" "sync" ) @@ -66,6 +67,12 @@ // function in order of the original heights values. func DoContours(tree *Tree, heights []float64, store func(*ContourResult)) { + if len(heights) == 0 { + return + } + + log.Printf("num heights: %d\n", len(heights)) + contours := make([]*ContourResult, len(heights)) for i, h := range heights {
--- a/pkg/octree/loader.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/loader.go Fri Mar 01 11:06:27 2019 +0100 @@ -111,7 +111,7 @@ return tree, nil } -func deserialize(data []byte) (*Tree, error) { +func Deserialize(data []byte) (*Tree, error) { return loadReader( bufio.NewReader( snappy.NewReader(
--- a/pkg/octree/tin.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/tin.go Fri Mar 01 11:06:27 2019 +0100 @@ -25,6 +25,7 @@ "math" "gemma.intevation.de/gemma/pkg/models" + "gemma.intevation.de/gemma/pkg/wkb" ) var ( @@ -62,9 +63,9 @@ switch { case err != nil: return err - case endian == wkbNDR: + case endian == wkb.NDR: order = binary.LittleEndian - case endian == wkbXDR: + case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) @@ -76,7 +77,7 @@ switch { case err != nil: return err - case geomType != wkbTinZ: + case geomType != wkb.TinZ: return fmt.Errorf("unknown geometry type %x", geomType) } @@ -113,9 +114,9 @@ switch { case err != nil: return err - case endian == wkbNDR: + case endian == wkb.NDR: order = binary.LittleEndian - case endian == wkbXDR: + case endian == wkb.XDR: order = binary.BigEndian default: return fmt.Errorf("unknown byte order %x", endian) @@ -125,7 +126,7 @@ switch { case err != nil: return err - case geomType != wkbTriangleZ: + case geomType != wkb.TriangleZ: return fmt.Errorf("unknown geometry type %d", geomType) }
--- a/pkg/octree/tree.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/tree.go Fri Mar 01 11:06:27 2019 +0100 @@ -32,6 +32,10 @@ Max Vertex } +func (ot *Tree) Vertices() []Vertex { + return ot.vertices +} + var scale = [4][4]float64{ {0.0, 0.0, 0.5, 0.5}, {0.5, 0.0, 1.0, 0.5}, @@ -39,6 +43,76 @@ {0.5, 0.5, 1.0, 1.0}, } +func (ot *Tree) Value(x, y float64) (float64, bool) { + + // out of bounding box + if x < ot.Min.X || ot.Max.X < x || + y < ot.Min.Y || ot.Max.Y < y { + 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}} + + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + + if top.pos > 0 { // node + if index := ot.index[top.pos:]; len(index) > 7 { + 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, + } + if nbox.Contains(x, y) { + if a != 0 { + stack = append(stack, frame{a, nbox}) + } + if b != 0 { + stack = append(stack, frame{b, nbox}) + } + break + } + } + } + } else { // leaf + pos := -top.pos - 1 + n := ot.index[pos] + indices := ot.index[pos+1 : pos+1+n] + + for _, idx := range indices { + tri := ot.triangles[idx] + t := Triangle{ + ot.vertices[tri[0]], + ot.vertices[tri[1]], + ot.vertices[tri[2]], + } + if t.Contains(x, y) { + return t.Plane3D().Z(x, y), true + } + } + } + } + + return 0, false +} + // Vertical does a vertical cross cut from (x1, y1) to (x2, y2). func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) { @@ -164,12 +238,14 @@ } else { max = mid } - if index := ot.index[pos:]; len(index) > 3 { - stack = append(stack, - frame{index[0], min, max}, - frame{index[1], min, max}, - frame{index[2], min, max}, - frame{index[3], min, max}) + if pos < int32(len(ot.index)) { + if index := ot.index[pos:]; len(index) > 3 { + stack = append(stack, + frame{index[0], min, max}, + frame{index[1], min, max}, + frame{index[2], min, max}, + frame{index[3], min, max}) + } } } else { // leaf pos = -pos - 1
--- a/pkg/octree/vertex.go Fri Mar 01 10:53:52 2019 +0100 +++ b/pkg/octree/vertex.go Fri Mar 01 11:06:27 2019 +0100 @@ -20,6 +20,8 @@ "log" "math" "sort" + + "gemma.intevation.de/gemma/pkg/wkb" ) type ( @@ -60,8 +62,54 @@ B float64 C float64 } + + Plane3D struct { + A float64 + B float64 + C float64 + D float64 + } ) +func (t *Triangle) Plane3D() Plane3D { + + v0 := t[1].Sub(t[0]) + v1 := t[2].Sub(t[0]) + n := v0.Cross(v1).Normalize() + + // x*nx+ y*ny+ z*nz + d = 0 + // d = - (x*nx+ y*ny + z*nz) + d := -t[0].Dot(n) + return Plane3D{ + A: n.X, + B: n.Y, + C: n.Z, + D: d, + } +} + +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 +} + +func (v Vertex) Normalize() Vertex { + s := 1 / v.Length() + return Vertex{ + X: s * v.X, + Y: s * v.Y, + Z: s * v.Z, + } +} + +func (v Vertex) Dot(w Vertex) float64 { + return v.X*w.X + v.Y*w.Y + v.Z*w.Z +} + +func (v Vertex) Length() float64 { + return math.Sqrt(v.Dot(v)) +} + // Minimize adjust this vertex v to hold the minimum // values component-wise of v and w. func (v *Vertex) Minimize(w Vertex) { @@ -166,6 +214,30 @@ return 0 } +func (v Vertex) Dot2(w Vertex) float64 { + return v.X*w.X + v.Y*w.Y +} + +func (t *Triangle) Contains(x, y float64) bool { + v0 := t[2].Sub(t[0]) + v1 := t[1].Sub(t[0]) + v2 := Vertex{X: x, Y: y}.Sub(t[0]) + + dot00 := v0.Dot2(v0) + dot01 := v0.Dot2(v1) + dot02 := v0.Dot2(v2) + dot11 := v1.Dot2(v1) + dot12 := v1.Dot2(v2) + + // Compute barycentric coordinates + invDenom := 1 / (dot00*dot11 - dot01*dot01) + u := (dot11*dot02 - dot01*dot12) * invDenom + v := (dot00*dot12 - dot01*dot02) * invDenom + + // Check if point is in triangle + return u >= 0 && v >= 0 && u+v < 1 +} + // IntersectHorizontal calculates the line string that // results when cutting a triangle a a certain height. // Can be empty (on intersection), @@ -351,13 +423,13 @@ buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbMultiLineStringZ) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.MultiLineStringZ) binary.Write(buf, binary.LittleEndian, uint32(len(mls))) for _, ml := range mls { - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbLineStringZ) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.LineStringZ) binary.Write(buf, binary.LittleEndian, uint32(len(ml))) for _, p := range ml { binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) @@ -381,13 +453,13 @@ buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbMultiLineString) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.MultiLineString) binary.Write(buf, binary.LittleEndian, uint32(len(mls))) for _, ml := range mls { - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbLineString) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.LineString) binary.Write(buf, binary.LittleEndian, uint32(len(ml))) for _, p := range ml { binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) @@ -518,6 +590,11 @@ a.Y2 < a.Y1 || a.Y2 < b.Y1) } +func (a Box2D) Contains(x, y float64) bool { + return a.X1 <= x && x <= a.X2 && + a.Y1 <= y && y <= a.Y2 +} + // Xi returns the i-th x component. func (a Box2D) Xi(i int) float64 { if i == 0 { @@ -843,13 +920,17 @@ buf := bytes.NewBuffer(make([]byte, 0, size)) - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbMultiPointZ) + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.MultiPointZ) binary.Write(buf, binary.LittleEndian, uint32(len(mpz))) + perPoint := bytes.NewBuffer(make([]byte, 0, 1+4)) + binary.Write(perPoint, binary.LittleEndian, wkb.NDR) + binary.Write(perPoint, binary.LittleEndian, wkb.PointZ) + hdr := perPoint.Bytes() + for _, p := range mpz { - binary.Write(buf, binary.LittleEndian, wkbNDR) - binary.Write(buf, binary.LittleEndian, wkbPointZ) + buf.Write(hdr) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z))
--- a/pkg/octree/wkb.go Fri Mar 01 10:53:52 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,28 +0,0 @@ -// 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 - -const ( - wkbXDR byte = 0 - wkbNDR byte = 1 - - wkbLineString uint32 = 2 - wkbMultiLineString uint32 = 5 - wkbPointZ uint32 = 1000 + 1 - wkbLineStringZ uint32 = 1000 + 2 - wkbMultiPointZ uint32 = 1000 + 4 - wkbMultiLineStringZ uint32 = 1000 + 5 - wkbTinZ uint32 = 1000 + 16 - wkbTriangleZ uint32 = 1000 + 17 -)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/wkb/wkb.go Fri Mar 01 11:06:27 2019 +0100 @@ -0,0 +1,30 @@ +// 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 wkb + +const ( + XDR byte = 0 + NDR byte = 1 + + Point uint32 = 1 + LineString uint32 = 2 + Polygon uint32 = 3 + MultiLineString uint32 = 5 + PointZ uint32 = 1000 + 1 + LineStringZ uint32 = 1000 + 2 + MultiPointZ uint32 = 1000 + 4 + MultiLineStringZ uint32 = 1000 + 5 + TinZ uint32 = 1000 + 16 + TriangleZ uint32 = 1000 + 17 +)