Mercurial > gemma
changeset 2493:e3487fa9284d octree-diff
Merged default into octree-diff branch.
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Mon, 04 Mar 2019 11:56:07 +0100 |
| parents | 10681749371d (diff) 408e0f4d4008 (current diff) |
| children | a727e0426240 |
| files | |
| diffstat | 20 files changed, 2127 insertions(+), 105 deletions(-) [+] |
line wrap: on
line diff
--- a/3rdpartylibs.sh Mon Mar 04 10:26:19 2019 +0100 +++ b/3rdpartylibs.sh Mon Mar 04 11:56:07 2019 +0100 @@ -35,6 +35,9 @@ go get -u -v github.com/robfig/cron # MIT +go get -u -v github.com/tidwall/rtree +# MIT + ## list of additional licenses that get fetched and installed as dependencies # github.com/fsnotify/fsnotify/ BSD-3-Clause # github.com/hashicorp/hcl/ MPL-2.0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/octreediff/db.go Mon Mar 04 11:56:07 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 Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,486 @@ +// 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 ( + "context" + "database/sql" + "errors" + "flag" + "fmt" + "log" + "runtime" + "sync" + "time" + + "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/octree" +) + +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` + + insertContourSQL = ` +INSERT INTO diff_contour_lines ( + height, + lines +) +SELECT + $4, + ST_Transform( + ST_Multi( + ST_CollectionExtract( + ST_SimplifyPreserveTopology( + ST_Multi(ST_Collectionextract( + ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 2)), + $3 + ), + 2 + ) + ), + 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 + 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 +) +INSERT INTO diff_contour_lines ( + height, + lines +) +SELECT + $5, + ST_Transform( + ST_Multi( + ST_intersection( + ST_CollectionExtract( + ST_SimplifyPreserveTopology( + ST_Multi(ST_Collectionextract( + ST_MakeValid(ST_GeomFromWKB($6, $2::integer)), 2)), + $7 + ), + 2 + ), + area + ) + ), + 4326 + ) + FROM inter` +) + +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) triangulate() (*octree.Triangulation, error) { + vertices := make([]octree.Vertex, len(pm)) + var i int + for p, z := range pm { + vertices[i] = octree.Vertex{X: p.x, Y: p.y, Z: z} + i++ + } + return octree.Triangulate(vertices) +} + +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 + } + defer tx.Rollback() + + 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)) + + 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 + } + + now = time.Now() + log.Printf("triangulation took %v\n", now.Sub(last)) + last = now + + builder := octree.NewBuilder(tri.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) + } + } + + log.Printf("num heights: %d\n", len(heights)) + + var dataSize int + + stmt, err := tx.PrepareContext(ctx, insertContourSQLClipped) + if err != nil { + return err + } + defer stmt.Close() + + log.Println("do contours") + octree.DoContours(tree, heights, func(res *octree.ContourResult) { + if err == nil && len(res.Lines) > 0 { + log.Printf("%.2f: lines: %d\n", res.Height, len(res.Lines)) + wkb := res.Lines.AsWKB2D() + dataSize += len(wkb) + _, err = stmt.ExecContext( + ctx, + bottleneck, + first.EPSG, + firstDate, + secondDate, + res.Height, + wkb, + 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 Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/imports/wkb.go Mon Mar 04 11:56:07 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 Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/models/cross.go Mon Mar 04 11:56:07 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 Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/builder.go Mon Mar 04 11:56:07 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 Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/cache.go Mon Mar 04 11:56:07 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 Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/contours.go Mon Mar 04 11:56:07 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 {
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/hull.go Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,81 @@ +// Copyright (C) 2018 Michael Fogleman +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +package octree + +import "sort" + +func cross2D(p, a, b Vertex) float64 { + return (a.X-p.X)*(b.Y-p.Y) - (a.Y-p.Y)*(b.X-p.X) +} + +// ConvexHull returns the convex hull of the provided points. +func ConvexHull(points []Vertex) []Vertex { + // copy points + pointsCopy := make([]Vertex, len(points)) + copy(pointsCopy, points) + points = pointsCopy + + // sort points + sort.Slice(points, func(i, j int) bool { + a := points[i] + b := points[j] + if a.X != b.X { + return a.X < b.X + } + return a.Y < b.Y + }) + + // filter nearly-duplicate points + distinctPoints := points[:0] + for i, p := range points { + if i > 0 && p.SquaredDistance2D(points[i-1]) < eps { + continue + } + distinctPoints = append(distinctPoints, p) + } + points = distinctPoints + + // find upper and lower portions + var U, L []Vertex + for _, p := range points { + for len(U) > 1 && cross2D(U[len(U)-2], U[len(U)-1], p) > 0 { + U = U[:len(U)-1] + } + for len(L) > 1 && cross2D(L[len(L)-2], L[len(L)-1], p) < 0 { + L = L[:len(L)-1] + } + U = append(U, p) + L = append(L, p) + } + + // reverse upper portion + for i, j := 0, len(U)-1; i < j; i, j = i+1, j-1 { + U[i], U[j] = U[j], U[i] + } + + // construct complete hull + if len(U) > 0 { + U = U[:len(U)-1] + } + if len(L) > 0 { + L = L[:len(L)-1] + } + return append(L, U...) +}
--- a/pkg/octree/loader.go Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/loader.go Mon Mar 04 11:56:07 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(
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/node.go Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,49 @@ +// Copyright (C) 2018 Michael Fogleman +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +package octree + +type node struct { + i int32 + t int32 + prev *node + next *node +} + +func newNode(nodes []node, i int32, prev *node) *node { + n := &nodes[i] + n.i = i + if prev == nil { + n.prev = n + n.next = n + } else { + n.next = prev.next + n.prev = prev + prev.next.prev = n + prev.next = n + } + return n +} + +func (n *node) remove() *node { + n.prev.next = n.next + n.next.prev = n.prev + n.i = -1 + return n.prev +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/polygon.go Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,345 @@ +// 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 + +import ( + "bytes" + "encoding/binary" + "fmt" + "math" + + "github.com/tidwall/rtree" + + "gemma.intevation.de/gemma/pkg/wkb" +) + +type ( + ring []float64 + + Polygon struct { + // TODO: Implement me! + rings []ring + + indices []*rtree.RTree + } + + IntersectionType byte + + lineSegment []float64 +) + +const ( + IntersectionInside IntersectionType = iota + IntersectionOutSide + IntersectionOverlaps +) + +func (ls lineSegment) Rect(interface{}) ([]float64, []float64) { + + var min, max [2]float64 + + if ls[0] < ls[2] { + min[0] = ls[0] + max[0] = ls[2] + } else { + min[0] = ls[2] + max[0] = ls[0] + } + + if ls[1] < ls[3] { + min[1] = ls[1] + max[1] = ls[3] + } else { + min[1] = ls[3] + max[1] = ls[1] + } + + return min[:], max[:] +} + +func (p *Polygon) Indexify() { + indices := make([]*rtree.RTree, len(p.rings)) + + for i := range indices { + index := rtree.New(nil) + indices[i] = index + + rng := p.rings[i] + + for i := 0; i < len(rng); i += 2 { + var ls lineSegment + if i+4 <= len(rng) { + ls = lineSegment(rng[i : i+4]) + } else { + ls = []float64{rng[i], rng[i+1], rng[0], rng[1]} + } + index.Insert(ls) + } + } + + p.indices = indices +} + +func (a Box2D) InterectsSegment(ls lineSegment) bool { + p1x := ls[0] + p1y := ls[1] + p2x := ls[2] + p2y := ls[3] + + left := a.X1 + right := a.X2 + top := a.Y1 + bottom := a.Y2 + + // The direction of the ray + dx := p2x - p1x + dy := p2y - p1y + + min, max := 0.0, 1.0 + + var t0, t1 float64 + + // Left and right sides. + // - If the line is parallel to the y axis. + if dx == 0 { + if p1x < left || p1x > right { + return false + } + } else { + // - Make sure t0 holds the smaller value by checking the direction of the line. + if dx > 0 { + t0 = (left - p1x) / dx + t1 = (right - p1x) / dx + } else { + t1 = (left - p1x) / dx + t0 = (right - p1x) / dx + } + + if t0 > min { + min = t0 + } + if t1 < max { + max = t1 + } + if min > max || max < 0 { + return false + } + } + + // The top and bottom side. + // - If the line is parallel to the x axis. + if dy == 0 { + if p1y < top || p1y > bottom { + return false + } + } else { + // - Make sure t0 holds the smaller value by checking the direction of the line. + if dy > 0 { + t0 = (top - p1y) / dy + t1 = (bottom - p1y) / dy + } else { + t1 = (top - p1y) / dy + t0 = (bottom - p1y) / dy + } + + if t0 > min { + min = t0 + } + if t1 < max { + max = t1 + } + if min > max || max < 0 { + return false + } + } + + // The point of intersection + // ix = p1x + dx*min + // iy = p1y + dy*min + return true +} + +func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType { + + if len(p.rings) == 0 { + return IntersectionOutSide + } + + for _, index := range p.indices { + var intersects bool + index.Search(box, func(item rtree.Item) bool { + ls := item.(lineSegment) + if box.InterectsSegment(ls) { + intersects = true + return false + } + return true + }) + if intersects { + return IntersectionOverlaps + } + } + + // No intersection -> check inside or outside + // if an abritrary point is inside or not. + point := []float64{box.X1, box.Y1} + + // Check holes first: inside a hole means outside. + if len(p.rings) > 0 { + for _, hole := range p.rings[1:] { + if hole.contains(point) { + return IntersectionOutSide + } + } + } + if p.rings[0].contains(point) { + return IntersectionInside + } + return IntersectionOutSide +} + +func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType { + // TODO: Implement me + return IntersectionOutSide +} + +func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 } + +func (rng ring) contains(point []float64) bool { + if !rng.isClosed() { + return false + } + + contains := intersectsWithRaycast(point, rng[:2], rng[len(rng)-2:len(rng)]) + + for i := 2; i < len(rng); i += 2 { + if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) { + contains = !contains + } + } + + return contains +} + +// Using the raycast algorithm, this returns whether or not the passed in point +// Intersects with the edge drawn by the passed in start and end points. +// Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go +func intersectsWithRaycast(point, start, end []float64) bool { + + // Always ensure that the the first point + // has a y coordinate that is less than the second point + if start[1] > end[1] { + // Switch the points if otherwise. + start, end = end, start + } + + // Move the point's y coordinate + // outside of the bounds of the testing region + // so we can start drawing a ray + for point[1] == start[1] || point[1] == end[1] { + y := math.Nextafter(point[1], math.Inf(1)) + point = []float64{point[0], y} + } + + // If we are outside of the polygon, indicate so. + if point[1] < start[1] || point[1] > end[1] { + return false + } + + if start[0] > end[0] { + if point[0] > start[0] { + return false + } + if point[0] < end[0] { + return true + } + } else { + if point[0] > end[0] { + return false + } + if point[0] < start[0] { + return true + } + } + + raySlope := (point[1] - start[1]) / (point[0] - start[0]) + diagSlope := (end[1] - start[1]) / (end[0] - start[0]) + + return raySlope >= diagSlope +} + +func (p *Polygon) FromWKB(data []byte) error { + + r := bytes.NewReader(data) + + endian, err := r.ReadByte() + + var order binary.ByteOrder + + switch { + case err != nil: + return err + case endian == wkb.NDR: + order = binary.LittleEndian + case endian == wkb.XDR: + order = binary.BigEndian + default: + return fmt.Errorf("unknown byte order %x", endian) + } + + var geomType uint32 + err = binary.Read(r, order, &geomType) + + switch { + case err != nil: + return err + case geomType != wkb.Polygon: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var numRings uint32 + if err = binary.Read(r, order, &numRings); err != nil { + return err + } + + rngs := make([]ring, numRings) + + for rng := uint32(0); rng < numRings; rng++ { + var numVertices uint32 + if err = binary.Read(r, order, &numVertices); err != nil { + return err + } + + numVertices *= 2 + vertices := make([]float64, numVertices) + + for v := uint32(0); v < numVertices; v += 2 { + var lat, lon uint64 + if err = binary.Read(r, order, &lat); err != nil { + return err + } + if err = binary.Read(r, order, &lon); err != nil { + return err + } + vertices[v] = math.Float64frombits(lat) + vertices[v+1] = math.Float64frombits(lon) + } + + rngs[rng] = vertices + } + + p.rings = rngs + + return nil +}
--- a/pkg/octree/tin.go Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/tin.go Mon Mar 04 11:56:07 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 Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/tree.go Mon Mar 04 11:56:07 2019 +0100 @@ -32,6 +32,15 @@ Max Vertex } +type boxFrame struct { + pos int32 + Box2D +} + +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 +48,170 @@ {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 + if x < ot.Min.X || ot.Max.X < x || + y < ot.Min.Y || ot.Max.Y < y { + return 0, false + } + + all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} + + stack := []boxFrame{{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, boxFrame{a, nbox}) + } + if b != 0 { + stack = append(stack, boxFrame{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)) { @@ -57,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] @@ -92,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}) } } } @@ -164,12 +332,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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/triangulation.go Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,116 @@ +// Copyright (C) 2018 Michael Fogleman +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +package octree + +import ( + "fmt" + "math" +) + +type Triangulation struct { + Points []Vertex + ConvexHull []Vertex + Triangles []int32 + Halfedges []int32 +} + +// Triangulate returns a Delaunay triangulation of the provided points. +func Triangulate(points []Vertex) (*Triangulation, error) { + t := newTriangulator(points) + err := t.triangulate() + return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err +} + +func (t *Triangulation) TriangleSlices() [][]int32 { + sl := make([][]int32, len(t.Triangles)/3) + var j int + for i := range sl { + sl[i] = t.Triangles[j : j+3] + j += 3 + } + return sl +} + +func (t *Triangulation) Tin() *Tin { + + min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} + max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} + + vertices := t.Points + + for _, v := range vertices { + min.Minimize(v) + max.Maximize(v) + } + + return &Tin{ + Vertices: vertices, + Triangles: t.TriangleSlices(), + Min: min, + Max: max, + } +} + +func (t *Triangulation) area() float64 { + var result float64 + points := t.Points + ts := t.Triangles + for i := 0; i < len(ts); i += 3 { + p0 := points[ts[i+0]] + p1 := points[ts[i+1]] + p2 := points[ts[i+2]] + result += area(p0, p1, p2) + } + return result / 2 +} + +// Validate performs several sanity checks on the Triangulation to check for +// potential errors. Returns nil if no issues were found. You normally +// shouldn't need to call this function but it can be useful for debugging. +func (t *Triangulation) Validate() error { + // verify halfedges + for i1, i2 := range t.Halfedges { + if i1 != -1 && t.Halfedges[i1] != i2 { + return fmt.Errorf("invalid halfedge connection") + } + if i2 != -1 && t.Halfedges[i2] != int32(i1) { + return fmt.Errorf("invalid halfedge connection") + } + } + + // verify convex hull area vs sum of triangle areas + hull1 := t.ConvexHull + hull2 := ConvexHull(t.Points) + area1 := polygonArea(hull1) + area2 := polygonArea(hull2) + area3 := t.area() + if math.Abs(area1-area2) > 1e-9 || math.Abs(area1-area3) > 1e-9 { + return fmt.Errorf("hull areas disagree: %f, %f, %f", area1, area2, area3) + } + + // verify convex hull perimeter + perimeter1 := polygonPerimeter(hull1) + perimeter2 := polygonPerimeter(hull2) + if math.Abs(perimeter1-perimeter2) > 1e-9 { + return fmt.Errorf("hull perimeters disagree: %f, %f", perimeter1, perimeter2) + } + + return nil +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/triangulator.go Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,375 @@ +// Copyright (C) 2018 Michael Fogleman +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +package octree + +import ( + "fmt" + "math" + "sort" +) + +type triangulator struct { + points []Vertex + squaredDistances []float64 + ids []int32 + center Vertex + triangles []int32 + halfedges []int32 + trianglesLen int + hull *node + hash []*node +} + +func newTriangulator(points []Vertex) *triangulator { + return &triangulator{points: points} +} + +// sorting a triangulator sorts the `ids` such that the referenced points +// are in order by their distance to `center` +func (a *triangulator) Len() int { + return len(a.points) +} + +func (a *triangulator) Swap(i, j int) { + a.ids[i], a.ids[j] = a.ids[j], a.ids[i] +} + +func (a *triangulator) Less(i, j int) bool { + d1 := a.squaredDistances[a.ids[i]] + d2 := a.squaredDistances[a.ids[j]] + if d1 != d2 { + return d1 < d2 + } + p1 := a.points[a.ids[i]] + p2 := a.points[a.ids[j]] + if p1.X != p2.X { + return p1.X < p2.X + } + return p1.Y < p2.Y +} + +func (tri *triangulator) triangulate() error { + points := tri.points + + n := len(points) + if n == 0 { + return nil + } + + tri.ids = make([]int32, n) + + // compute bounds + x0 := points[0].X + y0 := points[0].Y + x1 := points[0].X + y1 := points[0].Y + for i, p := range points { + if p.X < x0 { + x0 = p.X + } + if p.X > x1 { + x1 = p.X + } + if p.Y < y0 { + y0 = p.Y + } + if p.Y > y1 { + y1 = p.Y + } + tri.ids[i] = int32(i) + } + + var i0, i1, i2 int32 + + // pick a seed point close to midpoint + m := Vertex{X: (x0 + x1) / 2, Y: (y0 + y1) / 2} + minDist := infinity + for i, p := range points { + d := p.SquaredDistance2D(m) + if d < minDist { + i0 = int32(i) + minDist = d + } + } + + // find point closest to seed point + minDist = infinity + for i, p := range points { + if int32(i) == i0 { + continue + } + d := p.SquaredDistance2D(points[i0]) + if d > 0 && d < minDist { + i1 = int32(i) + minDist = d + } + } + + // find the third point which forms the smallest circumcircle + minRadius := infinity + for i, p := range points { + if int32(i) == i0 || int32(i) == i1 { + continue + } + r := circumradius(points[i0], points[i1], p) + if r < minRadius { + i2 = int32(i) + minRadius = r + } + } + if minRadius == infinity { + return fmt.Errorf("No Delaunay triangulation exists for this input.") + } + + // swap the order of the seed points for counter-clockwise orientation + if area(points[i0], points[i1], points[i2]) < 0 { + i1, i2 = i2, i1 + } + + tri.center = circumcenter(points[i0], points[i1], points[i2]) + + // sort the points by distance from the seed triangle circumcenter + tri.squaredDistances = make([]float64, n) + for i, p := range points { + tri.squaredDistances[i] = p.SquaredDistance2D(tri.center) + } + sort.Sort(tri) + + // initialize a hash table for storing edges of the advancing convex hull + hashSize := int(math.Ceil(math.Sqrt(float64(n)))) + tri.hash = make([]*node, hashSize) + + // initialize a circular doubly-linked list that will hold an advancing convex hull + nodes := make([]node, n) + + e := newNode(nodes, i0, nil) + e.t = 0 + tri.hashEdge(e) + + e = newNode(nodes, i1, e) + e.t = 1 + tri.hashEdge(e) + + e = newNode(nodes, i2, e) + e.t = 2 + tri.hashEdge(e) + + tri.hull = e + + maxTriangles := 2*n - 5 + tri.triangles = make([]int32, maxTriangles*3) + tri.halfedges = make([]int32, maxTriangles*3) + + tri.addTriangle(i0, i1, i2, -1, -1, -1) + + pp := Vertex{X: infinity, Y: infinity} + for k := 0; k < n; k++ { + i := tri.ids[k] + p := points[i] + + // skip nearly-duplicate points + if p.SquaredDistance2D(pp) < eps { + continue + } + pp = p + + // skip seed triangle points + if i == int32(i0) || i == int32(i1) || i == int32(i2) { + continue + } + + // find a visible edge on the convex hull using edge hash + var start *node + key := tri.hashKey(p) + for j := 0; j < len(tri.hash); j++ { + start = tri.hash[key] + if start != nil && start.i >= 0 { + break + } + key++ + if key >= len(tri.hash) { + key = 0 + } + } + start = start.prev + + e := start + for area(p, points[e.i], points[e.next.i]) >= 0 { + e = e.next + if e == start { + e = nil + break + } + } + if e == nil { + // likely a near-duplicate point; skip it + continue + } + walkBack := e == start + + // add the first triangle from the point + t := tri.addTriangle(e.i, i, e.next.i, -1, -1, e.t) + e.t = t // keep track of boundary triangles on the hull + e = newNode(nodes, i, e) + + // recursively flip triangles from the point until they satisfy the Delaunay condition + e.t = tri.legalize(t + 2) + + // walk forward through the hull, adding more triangles and flipping recursively + q := e.next + for area(p, points[q.i], points[q.next.i]) < 0 { + t = tri.addTriangle(q.i, i, q.next.i, q.prev.t, -1, q.t) + q.prev.t = tri.legalize(t + 2) + tri.hull = q.remove() + q = q.next + } + + if walkBack { + // walk backward from the other side, adding more triangles and flipping + q := e.prev + for area(p, points[q.prev.i], points[q.i]) < 0 { + t = tri.addTriangle(q.prev.i, i, q.i, -1, q.t, q.prev.t) + tri.legalize(t + 2) + q.prev.t = t + tri.hull = q.remove() + q = q.prev + } + } + + // save the two new edges in the hash table + tri.hashEdge(e) + tri.hashEdge(e.prev) + } + + tri.triangles = tri.triangles[:tri.trianglesLen] + tri.halfedges = tri.halfedges[:tri.trianglesLen] + + return nil +} + +func (t *triangulator) hashKey(point Vertex) int { + d := point.Sub2D(t.center) + return int(pseudoAngle(d.X, d.Y) * float64(len(t.hash))) +} + +func (t *triangulator) hashEdge(e *node) { + t.hash[t.hashKey(t.points[e.i])] = e +} + +func (t *triangulator) addTriangle(i0, i1, i2, a, b, c int32) int32 { + i := int32(t.trianglesLen) + t.triangles[i] = i0 + t.triangles[i+1] = i1 + t.triangles[i+2] = i2 + t.link(i, a) + t.link(i+1, b) + t.link(i+2, c) + t.trianglesLen += 3 + return i +} + +func (t *triangulator) link(a, b int32) { + t.halfedges[a] = b + if b >= 0 { + t.halfedges[b] = a + } +} + +func (t *triangulator) legalize(a int32) int32 { + // if the pair of triangles doesn't satisfy the Delaunay condition + // (p1 is inside the circumcircle of [p0, pl, pr]), flip them, + // then do the same check/flip recursively for the new pair of triangles + // + // pl pl + // /||\ / \ + // al/ || \bl al/ \a + // / || \ / \ + // / a||b \ flip /___ar___\ + // p0\ || /p1 => p0\---bl---/p1 + // \ || / \ / + // ar\ || /br b\ /br + // \||/ \ / + // pr pr + + b := t.halfedges[a] + + a0 := a - a%3 + b0 := b - b%3 + + al := a0 + (a+1)%3 + ar := a0 + (a+2)%3 + bl := b0 + (b+2)%3 + + if b < 0 { + return ar + } + + p0 := t.triangles[ar] + pr := t.triangles[a] + pl := t.triangles[al] + p1 := t.triangles[bl] + + illegal := inCircle(t.points[p0], t.points[pr], t.points[pl], t.points[p1]) + + if illegal { + t.triangles[a] = p1 + t.triangles[b] = p0 + + // edge swapped on the other side of the hull (rare) + // fix the halfedge reference + if t.halfedges[bl] == -1 { + e := t.hull + for { + if e.t == bl { + e.t = a + break + } + e = e.next + if e == t.hull { + break + } + } + } + + t.link(a, t.halfedges[bl]) + t.link(b, t.halfedges[ar]) + t.link(ar, bl) + + br := b0 + (b+1)%3 + + t.legalize(a) + return t.legalize(br) + } + + return ar +} + +func (t *triangulator) convexHull() []Vertex { + var result []Vertex + e := t.hull + for e != nil { + result = append(result, t.points[e.i]) + e = e.prev + if e == t.hull { + break + } + } + return result +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/util.go Mon Mar 04 11:56:07 2019 +0100 @@ -0,0 +1,37 @@ +// Copyright (C) 2018 Michael Fogleman +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +package octree + +import "math" + +var ( + eps = math.Nextafter(1, 2) - 1 + infinity = math.Inf(1) +) + +func pseudoAngle(dx, dy float64) float64 { + p := dx / (math.Abs(dx) + math.Abs(dy)) + if dy > 0 { + p = (3 - p) / 4 + } else { + p = (1 + p) / 4 + } + return math.Max(0, math.Min(1-eps, p)) +}
--- a/pkg/octree/vertex.go Mon Mar 04 10:26:19 2019 +0100 +++ b/pkg/octree/vertex.go Mon Mar 04 11:56:07 2019 +0100 @@ -20,6 +20,8 @@ "log" "math" "sort" + + "gemma.intevation.de/gemma/pkg/wkb" ) type ( @@ -60,8 +62,137 @@ 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)) +} + +func area(a, b, c Vertex) float64 { + return (b.Y-a.Y)*(c.X-b.X) - (b.X-a.X)*(c.Y-b.Y) +} + +func inCircle(a, b, c, p Vertex) bool { + dx := a.X - p.X + dy := a.Y - p.Y + ex := b.X - p.X + ey := b.Y - p.Y + fx := c.X - p.X + fy := c.Y - p.Y + + ap := dx*dx + dy*dy + bp := ex*ex + ey*ey + cp := fx*fx + fy*fy + + return dx*(ey*cp-bp*fy)-dy*(ex*cp-bp*fx)+ap*(ex*fy-ey*fx) < 0 +} + +func circumradius(a, b, c Vertex) float64 { + dx := b.X - a.X + dy := b.Y - a.Y + ex := c.X - a.X + ey := c.Y - a.Y + + bl := dx*dx + dy*dy + cl := ex*ex + ey*ey + d := dx*ey - dy*ex + + x := (ey*bl - dy*cl) * 0.5 / d + y := (dx*cl - ex*bl) * 0.5 / d + + r := x*x + y*y + + if bl == 0 || cl == 0 || d == 0 || r == 0 { + return infinity + } + + return r +} + +func circumcenter(a, b, c Vertex) Vertex { + dx := b.X - a.X + dy := b.Y - a.Y + ex := c.X - a.X + ey := c.Y - a.Y + + bl := dx*dx + dy*dy + cl := ex*ex + ey*ey + d := dx*ey - dy*ex + + x := a.X + (ey*bl-dy*cl)*0.5/d + y := a.Y + (dx*cl-ex*bl)*0.5/d + + return Vertex{X: x, Y: y} +} + +func polygonArea(points []Vertex) float64 { + var result float64 + for i, p := range points { + q := points[(i+1)%len(points)] + result += (p.X - q.X) * (p.Y + q.Y) + } + return result / 2 +} + +func polygonPerimeter(points []Vertex) float64 { + if len(points) == 0 { + return 0 + } + var result float64 + q := points[len(points)-1] + for _, p := range points { + result += p.Distance2D(q) + q = p + } + return result +} + +func (v Vertex) Distance2D(w Vertex) float64 { + return math.Hypot(v.X-w.X, v.Y-w.Y) +} + // Minimize adjust this vertex v to hold the minimum // values component-wise of v and w. func (v *Vertex) Minimize(w Vertex) { @@ -90,6 +221,20 @@ } } +func (v Vertex) SquaredDistance2D(w Vertex) float64 { + dx := v.X - w.X + dy := v.Y - w.Y + return dx*dx + dy*dy +} + +// Sub2D returns (v - w) component-wise. +func (v Vertex) Sub2D(w Vertex) Vertex { + return Vertex{ + X: v.X - w.X, + Y: v.Y - w.Y, + } +} + // Sub returns (v - w) component-wise. func (v Vertex) Sub(w Vertex) Vertex { return Vertex{ @@ -166,6 +311,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 +520,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 +550,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)) @@ -512,12 +681,21 @@ return out } +func (a Box2D) Rect(interface{}) ([]float64, []float64) { + return []float64{a.X1, a.Y1}, []float64{a.X2, a.Y2} +} + // Intersects checks if two Box2Ds intersect. func (a Box2D) Intersects(b Box2D) bool { return !(a.X2 < a.X1 || a.X2 < b.X1 || 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 +1021,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 Mon Mar 04 10:26:19 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 Mon Mar 04 11:56:07 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 +)