Mercurial > gemma
changeset 4827:f4abfd0ee8ad remove-octree-debris
Renamed octree package to mesh as there is no octree any more.
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Tue, 05 Nov 2019 14:30:22 +0100 |
| parents | ec5afada70ec |
| children | 39ee00d06309 |
| files | pkg/controllers/cross.go pkg/controllers/diff.go pkg/imports/isr.go pkg/imports/sr.go pkg/mesh/areas.go pkg/mesh/cache.go pkg/mesh/classbreaks.go pkg/mesh/hull.go pkg/mesh/loader.go pkg/mesh/node.go pkg/mesh/plane2d_test.go pkg/mesh/polygon.go pkg/mesh/raster.go pkg/mesh/strtree.go pkg/mesh/tin.go pkg/mesh/triangulation.go pkg/mesh/triangulator.go pkg/mesh/util.go pkg/mesh/vertex.go pkg/octree/areas.go pkg/octree/cache.go pkg/octree/classbreaks.go pkg/octree/hull.go pkg/octree/loader.go pkg/octree/node.go pkg/octree/plane2d_test.go pkg/octree/polygon.go pkg/octree/raster.go pkg/octree/strtree.go pkg/octree/tin.go pkg/octree/triangulation.go pkg/octree/triangulator.go pkg/octree/util.go pkg/octree/vertex.go |
| diffstat | 34 files changed, 4463 insertions(+), 4462 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/controllers/cross.go Tue Nov 05 13:01:37 2019 +0100 +++ b/pkg/controllers/cross.go Tue Nov 05 14:30:22 2019 +0100 @@ -21,8 +21,8 @@ "net/http" "time" + "gemma.intevation.de/gemma/pkg/mesh" "gemma.intevation.de/gemma/pkg/models" - "gemma.intevation.de/gemma/pkg/octree" mw "gemma.intevation.de/gemma/pkg/middleware" ) @@ -52,7 +52,7 @@ func projectBack( ctx context.Context, - line octree.MultiLineStringZ, + line mesh.MultiLineStringZ, epsg uint32, conn *sql.Conn, ) (models.GeoJSONMultiLineCoordinatesZ, error) { @@ -74,7 +74,7 @@ ctx := req.Context() conn := mw.JSONConn(req) - tree, err := octree.FromCache( + tree, err := mesh.FromCache( ctx, conn, csi.Properties.Bottleneck, csi.Properties.Date.Time) @@ -117,16 +117,16 @@ start = time.Now() - var segments octree.MultiLineStringZ + var segments mesh.MultiLineStringZ for i := 0; i < len(coords)-1; i++ { c1 := &coords[i] c2 := &coords[i+1] - verticalLine := octree.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) + verticalLine := mesh.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) - var line octree.MultiLineStringZ - tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *octree.Triangle) { + var line mesh.MultiLineStringZ + tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *mesh.Triangle) { if ls := verticalLine.Intersection(t); len(ls) > 0 { line = append(line, ls) }
--- a/pkg/controllers/diff.go Tue Nov 05 13:01:37 2019 +0100 +++ b/pkg/controllers/diff.go Tue Nov 05 14:30:22 2019 +0100 @@ -25,8 +25,8 @@ "gemma.intevation.de/gemma/pkg/auth" "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/mesh" "gemma.intevation.de/gemma/pkg/models" - "gemma.intevation.de/gemma/pkg/octree" mw "gemma.intevation.de/gemma/pkg/middleware" ) @@ -162,7 +162,7 @@ start := time.Now() begin := start - minuendTree, err := octree.FromCache( + minuendTree, err := mesh.FromCache( ctx, conn, dci.Bottleneck, dci.Minuend.Time) @@ -179,7 +179,7 @@ epsg := minuendTree.EPSG() - var box octree.Box2D + var box mesh.Box2D switch err := conn.QueryRowContext( ctx, @@ -203,7 +203,7 @@ box.X1, box.Y1, box.X2, box.Y2) start = time.Now() - raster := octree.NewRaster(box, isoCellSize) + raster := mesh.NewRaster(box, isoCellSize) raster.Rasterize(minuendTree.Value) log.Printf("info: rasterizing minuend took %v\n", time.Since(start)) @@ -211,7 +211,7 @@ start = time.Now() - subtrahendTree, err := octree.FromCache( + subtrahendTree, err := mesh.FromCache( ctx, conn, dci.Bottleneck, dci.Subtrahend.Time) @@ -253,15 +253,15 @@ var heights []float64 - heights, err = octree.LoadClassBreaks( + heights, err = mesh.LoadClassBreaks( ctx, tx, "morphology_classbreaks_compare") if err != nil { log.Printf("warn: Loading class breaks failed: %v\n", err) err = nil - heights = octree.SampleDiffHeights(zMin, zMax, contourStep) + heights = mesh.SampleDiffHeights(zMin, zMax, contourStep) } else { - heights = octree.ExtrapolateClassBreaks(heights, zMin, zMax) + heights = mesh.ExtrapolateClassBreaks(heights, zMin, zMax) } heights = common.DedupFloat64s(heights)
--- a/pkg/imports/isr.go Tue Nov 05 13:01:37 2019 +0100 +++ b/pkg/imports/isr.go Tue Nov 05 14:30:22 2019 +0100 @@ -19,7 +19,7 @@ "time" "gemma.intevation.de/gemma/pkg/common" - "gemma.intevation.de/gemma/pkg/octree" + "gemma.intevation.de/gemma/pkg/mesh" ) type IsoRefresh struct { @@ -123,7 +123,7 @@ time.Since(start)) }() - heights, err := octree.ParseClassBreaks(isr.ClassBreaks) + heights, err := mesh.ParseClassBreaks(isr.ClassBreaks) if err != nil { return nil, err } @@ -171,11 +171,11 @@ // For all sounding results in bottleneck. for _, sr := range bn.srs { - tree, err := octree.FetchMeshDirectly(ctx, tx, sr) + tree, err := mesh.FetchMeshDirectly(ctx, tx, sr) if err != nil { return err } - hs := octree.ExtrapolateClassBreaks(heights, tree.Min().Z, tree.Max().Z) + hs := mesh.ExtrapolateClassBreaks(heights, tree.Min().Z, tree.Max().Z) hs = common.DedupFloat64s(hs) // Delete the old iso areas. @@ -184,14 +184,14 @@ } // Calculate and store the iso areas. - box := octree.Box2D{ + box := mesh.Box2D{ X1: tree.Min().X, Y1: tree.Min().Y, X2: tree.Max().X, Y2: tree.Max().Y, } - raster := octree.NewRaster(box, isoCellSize) + raster := mesh.NewRaster(box, isoCellSize) raster.Rasterize(tree.Value) areas := raster.Trace(hs)
--- a/pkg/imports/sr.go Tue Nov 05 13:01:37 2019 +0100 +++ b/pkg/imports/sr.go Tue Nov 05 14:30:22 2019 +0100 @@ -35,8 +35,8 @@ shp "github.com/jonas-p/go-shp" "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/mesh" "gemma.intevation.de/gemma/pkg/models" - "gemma.intevation.de/gemma/pkg/octree" "gemma.intevation.de/gemma/pkg/wkb" ) @@ -314,8 +314,8 @@ ldc /= 100 xform = chainTransforms( xform, - func(v octree.Vertex) octree.Vertex { - return octree.Vertex{X: v.X, Y: v.Y, Z: v.Z + ldc} + func(v mesh.Vertex) mesh.Vertex { + return mesh.Vertex{X: v.X, Y: v.Y, Z: v.Z + ldc} }) m.DepthReference = depthReference } @@ -324,7 +324,7 @@ return nil, common.ToError(err) } - var xyz octree.MultiPointZ + var xyz mesh.MultiPointZ if z != nil { // Scanning ZIP file for *.xyz file. var xyzf *zip.File @@ -392,7 +392,7 @@ feedback Feedback, importID int64, m *models.SoundingResultMeta, - xyz octree.MultiPointZ, + xyz mesh.MultiPointZ, boundary polygonSlice, ) (interface{}, error) { @@ -429,7 +429,7 @@ feedback.Info("Triangulate XYZ data.") start = time.Now() - tri, err := octree.Triangulate(xyz) + tri, err := mesh.Triangulate(xyz) if err != nil { return nil, err } @@ -437,12 +437,12 @@ feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) var ( - clippingPolygon octree.Polygon - clippingPolygonBuffered octree.Polygon + clippingPolygon mesh.Polygon + clippingPolygonBuffered mesh.Polygon removed map[int32]struct{} polygonArea float64 clippingPolygonWKB []byte - tin *octree.Tin + tin *mesh.Tin ) if boundary == nil { @@ -450,7 +450,7 @@ tooLongEdge := tri.EstimateTooLong() feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) - var polygon octree.LineStringZ + var polygon mesh.LineStringZ start = time.Now() polygon, removed = tri.ConcaveHull(tooLongEdge) @@ -503,7 +503,7 @@ tin := tri.Tin() tin.EPSG = epsg - var str octree.STRTree + var str mesh.STRTree str.Build(tin) removed = str.Clip(&clippingPolygon) @@ -529,7 +529,7 @@ start = time.Now() tin := tri.Tin() - var virtual octree.STRTree + var virtual mesh.STRTree virtual.BuildWithout(tin, removed) feedback.Info("Building took %v", time.Since(start)) @@ -541,13 +541,13 @@ start = time.Now() - generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) + generated := make(mesh.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) - octree.GenerateRandomVertices( + mesh.GenerateRandomVertices( numPoints, tin.Min, tin.Max, virtual.Value, - func(vertices []octree.Vertex) { + func(vertices []mesh.Vertex) { generated = append(generated, vertices...) }) @@ -562,15 +562,15 @@ } dupes[key] = struct{}{} if z, ok := virtual.Value(x, y); ok { - generated = append(generated, octree.Vertex{X: x, Y: y, Z: z}) + generated = append(generated, mesh.Vertex{X: x, Y: y, Z: z}) } }) feedback.Info("Triangulate new point cloud.") - xyz = octree.MultiPointZ(generated) + xyz = mesh.MultiPointZ(generated) start = time.Now() - tri, err = octree.Triangulate(xyz) + tri, err = mesh.Triangulate(xyz) if err != nil { return nil, err } @@ -585,7 +585,7 @@ tin = tri.Tin() tin.EPSG = epsg - var str octree.STRTree + var str mesh.STRTree str.Build(tin) feedback.Info("Building clipping index took %v", time.Since(start)) @@ -597,7 +597,7 @@ feedback.Info("Number of triangles to clip %d.", len(removed)) start = time.Now() - final := octree.STRTree{Entries: 16} + final := mesh.STRTree{Entries: 16} final.BuildWithout(tin, removed) feedback.Info("Building final mesh took %v.", time.Since(start)) @@ -742,20 +742,20 @@ return &m, nil } -type vertexTransform func(octree.Vertex) octree.Vertex +type vertexTransform func(mesh.Vertex) mesh.Vertex -func identityTransform(v octree.Vertex) octree.Vertex { return v } +func identityTransform(v mesh.Vertex) mesh.Vertex { return v } -func negateZTransform(v octree.Vertex) octree.Vertex { - return octree.Vertex{X: v.X, Y: v.Y, Z: -v.Z} +func negateZTransform(v mesh.Vertex) mesh.Vertex { + return mesh.Vertex{X: v.X, Y: v.Y, Z: -v.Z} } func chainTransforms(a, b vertexTransform) vertexTransform { - return func(v octree.Vertex) octree.Vertex { return b(a(v)) } + return func(v mesh.Vertex) mesh.Vertex { return b(a(v)) } } -func loadXYZReader(r io.Reader, feedback Feedback, xform vertexTransform) (octree.MultiPointZ, error) { - mpz := make(octree.MultiPointZ, 0, 250000) +func loadXYZReader(r io.Reader, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { + mpz := make(mesh.MultiPointZ, 0, 250000) s := bufio.NewScanner(r) warnLimiter := common.WarningLimiter{Log: feedback.Warn, MaxWarnings: 100} @@ -764,7 +764,7 @@ for line := 1; s.Scan(); line++ { text := s.Text() - var p octree.Vertex + var p mesh.Vertex // fmt.Sscanf(text, "%f,%f,%f") is 4 times slower. idx := strings.IndexByte(text, ',') if idx == -1 { @@ -801,7 +801,7 @@ return mpz, nil } -func loadXYZ(f *zip.File, feedback Feedback, xform vertexTransform) (octree.MultiPointZ, error) { +func loadXYZ(f *zip.File, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { r, err := f.Open() if err != nil { return nil, err @@ -810,7 +810,7 @@ return loadXYZReader(r, feedback, xform) } -func loadXYZFile(f string, feedback Feedback, xform vertexTransform) (octree.MultiPointZ, error) { +func loadXYZFile(f string, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { r, err := os.Open(f) if err != nil { return nil, err @@ -859,11 +859,11 @@ ctx context.Context, tx *sql.Tx, feedback Feedback, - tree *octree.STRTree, + tree *mesh.STRTree, id int64, ) error { - heights, err := octree.LoadClassBreaks( + heights, err := mesh.LoadClassBreaks( ctx, tx, "morphology_classbreaks", ) @@ -879,7 +879,7 @@ heights = append(heights, h) } } else { - heights = octree.ExtrapolateClassBreaks(heights, minZ, maxZ) + heights = mesh.ExtrapolateClassBreaks(heights, minZ, maxZ) } /* @@ -898,7 +898,7 @@ ctx context.Context, tx *sql.Tx, feedback Feedback, - tree *octree.STRTree, + tree *mesh.STRTree, heights []float64, id int64, ) error { @@ -909,14 +909,14 @@ time.Since(total)) }() - box := octree.Box2D{ + box := mesh.Box2D{ X1: tree.Min().X, Y1: tree.Min().Y, X2: tree.Max().X, Y2: tree.Max().Y, } - raster := octree.NewRaster(box, isoCellSize) + raster := mesh.NewRaster(box, isoCellSize) raster.Rasterize(tree.Value) areas := raster.Trace(heights)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/areas.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,88 @@ +// 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 mesh + +import ( + "runtime" + "sync" + + "gemma.intevation.de/gemma/pkg/common" +) + +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 +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/cache.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,199 @@ +// 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 mesh + +import ( + "context" + "database/sql" + "sync" + "time" +) + +type ( + cacheKey struct { + date time.Time + bottleneck string + } + + cacheEntry struct { + checksum string + tree *STRTree + access time.Time + } + + // Cache holds Octrees for a defined amount of time in memory + // before they are released. + Cache struct { + sync.Mutex + entries map[cacheKey]*cacheEntry + } +) + +const ( + cleanupCacheSleep = 6 * time.Minute + maxCacheAge = 5 * time.Minute + maxCacheEntries = 4 +) + +const ( + directMeshSQL = ` +SELECT mesh_index FROM waterway.sounding_results +WHERE id = $1 +` + fetchMeshSQL = ` +SELECT mesh_checksum, mesh_index +FROM waterway.sounding_results +WHERE bottleneck_id = $1 AND date_info = $2::date + AND mesh_checksum IS NOT NULL AND mesh_index IS NOT NULL +` + checkMeshSQL = ` +SELECT CASE + WHEN mesh_checksum = $3 THEN NULL + ELSE mesh_index + END +FROM waterway.sounding_results +WHERE bottleneck_id = $1 AND date_info = $2::date + AND mesh_checksum IS NOT NULL AND mesh_index IS NOT NULL +` +) + +var cache = Cache{ + entries: map[cacheKey]*cacheEntry{}, +} + +func init() { + go cache.background() +} + +func (c *Cache) background() { + for { + time.Sleep(cleanupCacheSleep) + c.cleanup() + } +} + +func (c *Cache) cleanup() { + c.Lock() + defer c.Unlock() + good := time.Now().Add(-maxCacheAge) + for k, v := range c.entries { + if v.access.Before(good) { + delete(c.entries, k) + } + } +} + +// FromCache fetches an Octree from the global Octree cache. +func FromCache( + ctx context.Context, + conn *sql.Conn, + bottleneck string, date time.Time, +) (*STRTree, error) { + return cache.get(ctx, conn, bottleneck, date) +} + +// FetchOctreeDirectly loads an octree directly from the database. +func FetchMeshDirectly( + ctx context.Context, + tx *sql.Tx, + id int64, +) (*STRTree, error) { + var data []byte + err := tx.QueryRowContext(ctx, directMeshSQL, id).Scan(&data) + if err != nil { + return nil, err + } + 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, +) (*STRTree, error) { + c.Lock() + defer c.Unlock() + + key := cacheKey{date, bottleneck} + entry := c.entries[key] + + var data []byte + var checksum string + + if entry == nil { + // fetch from database + err := conn.QueryRowContext( + ctx, fetchMeshSQL, bottleneck, date).Scan(&checksum, &data) + switch { + case err == sql.ErrNoRows: + return nil, nil + case err != nil: + return nil, err + } + } else { + // check if we are not outdated. + err := conn.QueryRowContext( + ctx, checkMeshSQL, bottleneck, date, entry.checksum).Scan(&data) + switch { + case err == sql.ErrNoRows: + return nil, nil + case err != nil: + return nil, err + } + if data == nil { // we are still current + entry.access = time.Now() + return entry.tree, nil + } + } + + tree := new(STRTree) + + if err := tree.FromBytes(data); err != nil { + return nil, err + } + + now := time.Now() + + if entry != nil { + entry.tree = tree + entry.access = now + return tree, nil + } + + for len(c.entries) >= maxCacheEntries { + // Evict the entry that is accessed the longest time ago. + var oldestKey cacheKey + oldest := now + + for k, v := range c.entries { + if v.access.Before(oldest) { + oldest = v.access + oldestKey = k + } + } + delete(c.entries, oldestKey) + } + + c.entries[key] = &cacheEntry{ + checksum: checksum, + tree: tree, + access: now, + } + + return tree, nil +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/classbreaks.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,150 @@ +// 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 mesh + +import ( + "context" + "database/sql" + "errors" + "math" + "sort" + "strconv" + "strings" +) + +const ( + selectClassBreaksSQL = ` +SELECT config_val FROM sys_admin.system_config +WHERE config_key = $1` +) + +func SampleDiffHeights(min, max, step float64) []float64 { + var heights []float64 + switch { + case min >= 0: // All values positive. + for v := 0.0; v <= max; v += step { + if v >= min { + heights = append(heights, v) + } + } + case max <= 0: // All values negative. + for v := 0.0; v >= min; v -= step { + if v <= max { + heights = append(heights, v) + } + } + default: // Positive and negative. + for v := step; v <= max; v += step { + 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 >= min; v -= step { + heights = append(heights, v) + } + } + return heights +} + +func ParseClassBreaks(config string) ([]float64, error) { + + parts := strings.Split(config, ",") + classes := make([]float64, 0, len(parts)) + for _, part := range parts { + if idx := strings.IndexRune(part, ':'); idx >= 0 { + part = part[:idx] + } + if part = strings.TrimSpace(part); part == "" { + continue + } + v, err := strconv.ParseFloat(part, 64) + if err != nil { + return nil, err + } + classes = append(classes, v) + } + + sort.Float64s(classes) + return classes, nil +} + +func LoadClassBreaks(ctx context.Context, tx *sql.Tx, key string) ([]float64, error) { + + var config sql.NullString + + err := tx.QueryRowContext(ctx, selectClassBreaksSQL, key).Scan(&config) + + switch { + case err == sql.ErrNoRows: + return nil, nil + case err != nil: + return nil, err + case !config.Valid: + return nil, errors.New("Invalid config string") + } + + return ParseClassBreaks(config.String) +} + +func round(v float64) float64 { + return math.Round(v*10000) / 10000 +} + +func ExtrapolateClassBreaks(cbs []float64, min, max float64) []float64 { + if min > max { + min, max = max, min + } + + n := make([]float64, len(cbs)) + copy(n, cbs) + sort.Float64s(n) + + for len(n) > 0 && n[0] < min { + n = n[1:] + } + + if len(n) == 0 { + return n + } + + for len(n) > 0 && n[len(n)-1] > max { + n = n[:len(n)-1] + } + + if len(n) == 0 { + return n + } + + for min < n[0] { + diff := n[1] - n[0] + if diff == 0 { + break + } + m := make([]float64, len(n)+1) + m[0] = round(n[0] - diff) + copy(m[1:], n) + n = m + } + + for max > n[len(n)-1] { + diff := n[len(n)-1] - n[len(n)-2] + if diff == 0 { + break + } + n = append(n, round(n[len(n)-1]+diff)) + } + + return n +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/hull.go Tue Nov 05 14:30:22 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 mesh + +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...) +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/loader.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,166 @@ +// 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 mesh + +import ( + "bufio" + "bytes" + "compress/gzip" + "encoding/binary" + "log" +) + +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 + } + + 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) + } + } + + for i := range bboxes { + read(&bboxes[i].X1) + read(&bboxes[i].Y1) + read(&bboxes[i].X2) + read(&bboxes[i].Y2) + } + + return 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 + } + + 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", + 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 err + } + + log.Printf("info: vertices: %d\n", numVertices) + + vertices := make([]Vertex, numVertices) + t.Vertices = vertices + + for i := range vertices { + if err := vertices[i].Read(r); err != nil { + return err + } + } + + var numTriangles uint32 + if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil { + return err + } + + log.Printf("info: triangles: %d\n", numTriangles) + + indices := make([]int32, 3*numTriangles) + triangles := make([][]int32, numTriangles) + t.Triangles = triangles + + var last int32 + + for i := range triangles { + tri := indices[:3] + indices = indices[3:] + triangles[i] = tri + for j := range tri { + v, err := binary.ReadVarint(r) + if err != nil { + return err + } + value := int32(v) + last + tri[j] = value + last = value + } + } + + return nil +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/node.go Tue Nov 05 14:30:22 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 mesh + +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/mesh/plane2d_test.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,47 @@ +// 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 mesh + +import ( + "math" + "testing" +) + +func TestIntersection(t *testing.T) { + + table := []struct { + a [4]float64 + b [4]float64 + intersects bool + x, y float64 + }{ + {[4]float64{-1, -1, 1, 1}, [4]float64{-1, 1, 1, -1}, true, 0, 0}, + {[4]float64{0, 0, 1, 1}, [4]float64{0, 1, 1, 0}, true, 0.5, 0.5}, + {[4]float64{0, 0, 1, 0}, [4]float64{0, 1, 1, 1}, false, 0, 0}, + } + + for _, e := range table { + p1 := NewPlane2D(e.a[0], e.a[1], e.a[2], e.a[3]) + p2 := NewPlane2D(e.b[0], e.b[1], e.b[2], e.b[3]) + x, y, intersects := p1.Intersection(p2) + if intersects != e.intersects { + t.Fatalf("Have %t want %t\n", intersects, e.intersects) + } + if e.intersects { + if math.Abs(e.x-x) > epsPlane || math.Abs(e.y-y) > epsPlane { + t.Fatalf("Have (%f, %f)t want (%f, %f)\n", x, y, e.x, e.y) + } + } + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/polygon.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,505 @@ +// 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 mesh + +import ( + "bytes" + "encoding/binary" + "fmt" + "log" + "math" + + "github.com/tidwall/rtree" + + "gemma.intevation.de/gemma/pkg/wkb" +) + +type ( + ring []float64 + + Polygon struct { + 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 (ls lineSegment) intersects(a Box2D) 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 (ls lineSegment) intersectsLineSegment(o lineSegment) bool { + + p0 := ls[:2] + p1 := ls[2:4] + p2 := o[:2] + p3 := o[2:4] + + s10x := p1[0] - p0[0] + s10y := p1[1] - p0[1] + s32x := p3[0] - p2[0] + s32y := p3[1] - p2[1] + + den := s10x*s32y - s32x*s10y + + if den == 0 { + return false + } + + denPos := den > 0 + + s02x := p0[0] - p2[0] + s02y := p0[1] - p2[1] + + sNum := s10x*s02y - s10y*s02x + if sNum < 0 == denPos { + return false + } + + tNum := s32x*s02y - s32y*s02x + if tNum < 0 == denPos { + return false + } + + if sNum > den == denPos || tNum > den == denPos { + return false + } + + // t := tNum / den + // intersection at( p0[0] + (t * s10x), p0[1] + (t * s10y) ) + 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 { + if item.(lineSegment).intersects(box) { + intersects = true + return false + } + return true + }) + if intersects { + return IntersectionOverlaps + } + } + + // No intersection -> check inside or outside + // if an abritrary point is inside or not. + + // Check holes first: inside a hole means outside. + if len(p.rings) > 1 { + for _, hole := range p.rings[1:] { + if contains(hole, box.X1, box.Y1) { + return IntersectionOutSide + } + } + } + + // Check shell + if contains(p.rings[0], box.X1, box.Y1) { + return IntersectionInside + } + return IntersectionOutSide +} + +func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType { + box := t.BBox() + for _, index := range p.indices { + var intersects bool + index.Search(box, func(item rtree.Item) bool { + ls := item.(lineSegment) + other := make(lineSegment, 4) + for i := range t { + other[0] = t[i].X + other[1] = t[i].Y + other[2] = t[(i+1)%len(t)].X + other[3] = t[(i+1)%len(t)].Y + if ls.intersectsLineSegment(other) { + intersects = true + return false + } + } + return true + }) + if intersects { + return IntersectionOverlaps + } + } + // No intersection -> check inside or outside + // if an abritrary point is inside or not. + pX, pY := t[0].X, t[0].Y + + // Check holes first: inside a hole means outside. + if len(p.rings) > 1 { + for _, hole := range p.rings[1:] { + if contains(hole, pX, pY) { + return IntersectionOutSide + } + } + } + + // Check shell + if contains(p.rings[0], pX, pY) { + return IntersectionInside + } + return IntersectionOutSide +} + +func (rng ring) length() int { + return len(rng) / 2 +} + +func (rng ring) point(i int) (float64, float64) { + i *= 2 + return rng[i], rng[i+1] +} + +type segments interface { + length() int + point(int) (float64, float64) +} + +func contains(s segments, pX, pY float64) bool { + + n := s.length() + if n < 3 { + return false + } + + sX, sY := s.point(0) + eX, eY := s.point(n - 1) + + const eps = 0.0000001 + + if math.Abs(sX-eX) > eps || math.Abs(sY-eY) > eps { + // It's not closed! + return false + } + + var inside bool + + for i := 1; i < n; i++ { + eX, eY := s.point(i) + if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) { + inside = !inside + } + sX, sY = eX, eY + } + + return inside +} + +// 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(pX, pY, sX, sY, eX, eY float64) bool { + + // Always ensure that the the first point + // has a y coordinate that is less than the second point + if sY > eY { + // Switch the points if otherwise. + sX, sY, eX, eY = eX, eY, sX, sY + } + + // Move the point's y coordinate + // outside of the bounds of the testing region + // so we can start drawing a ray + for pY == sY || pY == eY { + pY = math.Nextafter(pY, math.Inf(1)) + } + + // If we are outside of the polygon, indicate so. + if pY < sY || pY > eY { + return false + } + + if sX > eX { + if pX > sX { + return false + } + if pX < eX { + return true + } + } else { + if pX > eX { + return false + } + if pX < sX { + return true + } + } + + raySlope := (pY - sY) / (pX - sX) + diagSlope := (eY - sY) / (eX - sX) + + return raySlope >= diagSlope +} + +func (p *Polygon) NumVertices(ring int) int { + if ring < 0 || ring >= len(p.rings) { + return 0 + } + return len(p.rings[ring]) / 2 +} + +func (p *Polygon) Vertices(ring int, fn func(float64, float64)) { + if ring < 0 || ring >= len(p.rings) { + return + } + rng := p.rings[ring] + for i := 0; i < len(rng); i += 2 { + fn(rng[i+0], rng[i+1]) + } +} + +func (p *Polygon) AsWKB() []byte { + + size := 1 + 4 + 4 + for _, r := range p.rings { + size += 4 + len(r)*8 + } + + buf := bytes.NewBuffer(make([]byte, 0, size)) + + binary.Write(buf, binary.LittleEndian, wkb.NDR) + binary.Write(buf, binary.LittleEndian, wkb.Polygon) + binary.Write(buf, binary.LittleEndian, uint32(len(p.rings))) + + for _, r := range p.rings { + binary.Write(buf, binary.LittleEndian, uint32(len(r)/2)) + for i := 0; i < len(r); i += 2 { + binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0])) + binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1])) + } + } + + return buf.Bytes() +} + +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) + + log.Printf("info: Number of rings: %d\n", len(rngs)) + + for rng := uint32(0); rng < numRings; rng++ { + var numVertices uint32 + if err = binary.Read(r, order, &numVertices); err != nil { + return err + } + + log.Printf("info: Number of vertices in ring %d: %d\n", rng, numVertices) + + 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 +} + +func (p *Polygon) FromLineStringZ(ls LineStringZ) { + r := make([]float64, 2*len(ls)) + var pos int + for i := range ls { + r[pos+0] = ls[i].X + r[pos+1] = ls[i].Y + pos += 2 + } + p.rings = []ring{r} +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/raster.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,405 @@ +// 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 mesh + +import ( + "log" + "math" + "runtime" + "sync" + "time" + + "github.com/fogleman/contourmap" + + "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/wkb" +) + +type Raster struct { + BBox Box2D + CellSize float64 + XCells int + YCells int + Cells []float64 +} + +const noData = -math.MaxFloat64 + +func NewRaster(bbox Box2D, cellSize float64) *Raster { + + width, height := bbox.Size() + + log.Printf("info raster extent: %.2f / %.2f", width, height) + + xCells := int(math.Ceil(width / cellSize)) + yCells := int(math.Ceil(height / cellSize)) + + log.Printf("info raster size: %d / %d\n", xCells, yCells) + + size := (xCells + 2) * (yCells + 2) + cells := make([]float64, size) + for i := range cells { + cells[i] = noData + } + return &Raster{ + BBox: bbox, + CellSize: cellSize, + XCells: xCells, + YCells: yCells, + Cells: cells, + } +} + +func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) { + var wg sync.WaitGroup + + rows := make(chan int) + + rasterRow := func() { + defer wg.Done() + quat := 0.25 * r.CellSize + for i := range rows { + pos := (i+1)*(r.XCells+2) + 1 + row := r.Cells[pos : pos+r.XCells] + py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 + px := r.BBox.X1 + r.CellSize/2 + y1 := py - quat + y2 := py + quat + for j := range row { + var n int + var sum float64 + + if v, ok := eval(px-quat, y1); ok { + sum = v + n = 1 + } + if v, ok := eval(px-quat, y2); ok { + sum += v + n++ + } + if v, ok := eval(px+quat, y1); ok { + sum += v + n++ + } + if v, ok := eval(px+quat, y2); ok { + sum += v + n++ + } + + if n > 0 { + row[j] = sum / float64(n) + } + px += r.CellSize + } + } + } + + for n := runtime.NumCPU(); n >= 1; n-- { + wg.Add(1) + go rasterRow() + } + + for i := 0; i < r.YCells; i++ { + rows <- i + } + close(rows) +} + +func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) { + var wg sync.WaitGroup + + rows := make(chan int) + + rasterRow := func() { + defer wg.Done() + quat := 0.25 * r.CellSize + for i := range rows { + pos := (i+1)*(r.XCells+2) + 1 + row := r.Cells[pos : pos+r.XCells] + py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 + px := r.BBox.X1 + r.CellSize/2 + y1 := py - quat + y2 := py + quat + for j, old := range row { + // only diff where values + if old == noData { + px += r.CellSize + continue + } + var n int + var sum float64 + + if v, ok := eval(px-quat, y1); ok { + sum = v + n = 1 + } + if v, ok := eval(px-quat, y2); ok { + sum += v + n++ + } + if v, ok := eval(px+quat, y1); ok { + sum += v + n++ + } + if v, ok := eval(px+quat, y2); ok { + sum += v + n++ + } + + if n > 0 { + row[j] -= sum / float64(n) + } else { + row[j] = noData + } + + px += r.CellSize + } + } + } + + for n := runtime.NumCPU(); n >= 1; n-- { + wg.Add(1) + go rasterRow() + } + + for i := 0; i < r.YCells; i++ { + rows <- i + } + close(rows) +} + +func (r *Raster) ZExtent() (float64, float64, bool) { + min, max := math.MaxFloat64, -math.MaxFloat64 + for _, v := range r.Cells { + if v == noData { + continue + } + if v < min { + min = v + } + if v > max { + max = v + } + } + return min, max, min != math.MaxFloat64 +} + +func (r *Raster) Trace(heights []float64) []wkb.MultiPolygonGeom { + start := time.Now() + + tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells) + + areas := make([]wkb.MultiPolygonGeom, len(heights)) + + reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize) + reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize) + + var wg sync.WaitGroup + + cnts := make(chan int) + + doContours := func() { + defer wg.Done() + for hIdx := range cnts { + if c := tracer.Contours(heights[hIdx]); len(c) > 0 { + areas[hIdx] = buildMultipolygon(c, reprojX, reprojY) + } + } + } + + for n := runtime.NumCPU(); n >= 1; n-- { + wg.Add(1) + go doContours() + } + + for i := range heights { + cnts <- i + } + close(cnts) + + wg.Wait() + log.Printf("info: Tracing areas took %v\n", time.Since(start)) + + return areas +} + +type contour []contourmap.Point + +type bboxNode struct { + box Box2D + cnt contour + children []*bboxNode +} + +func (cnt contour) contains(o contour) bool { + return contains(cnt, o[0].X, o[0].Y) || + contains(cnt, o[len(o)/2].X, o[len(o)/2].Y) +} + +func (cnt contour) length() int { + return len(cnt) +} + +func (cnt contour) point(i int) (float64, float64) { + return cnt[i].X, cnt[i].Y +} + +func (cnt contour) bbox() Box2D { + minX, minY := math.MaxFloat64, math.MaxFloat64 + maxX, maxY := -math.MaxFloat64, -math.MaxFloat64 + + for _, p := range cnt { + if p.X < minX { + minX = p.X + } + if p.X > maxX { + maxX = p.X + } + if p.Y < minY { + minY = p.Y + } + if p.Y > maxY { + maxY = p.Y + } + } + return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY} +} + +func (bn *bboxNode) insert(cnt contour, box Box2D) { + // check if children are inside new + var nr *bboxNode + + for i, r := range bn.children { + if r.box.Inside(box) && cnt.contains(r.cnt) { + if nr == nil { + nr = &bboxNode{box: box, cnt: cnt} + } + nr.children = append(nr.children, r) + bn.children[i] = nil + } + } + + // we have a new child + if nr != nil { + // compact the list + for i := len(bn.children) - 1; i >= 0; i-- { + if bn.children[i] == nil { + if i < len(bn.children)-1 { + copy(bn.children[i:], bn.children[i+1:]) + } + bn.children[len(bn.children)-1] = nil + bn.children = bn.children[:len(bn.children)-1] + } + } + bn.children = append(bn.children, nr) + return + } + + // check if new is inside an old + for _, r := range bn.children { + if box.Inside(r.box) && r.cnt.contains(cnt) { + r.insert(cnt, box) + return + } + } + + // its a new child node. + nr = &bboxNode{box: box, cnt: cnt} + bn.children = append(bn.children, nr) +} + +func (bn *bboxNode) insertRoot(cnt contour) { + bn.insert(cnt, cnt.bbox()) +} + +type bboxOutFunc func(contour, []contour) + +func (bn *bboxNode) generate(out bboxOutFunc) { + + var grands []*bboxNode + + holes := make([]contour, len(bn.children)) + + for i, ch := range bn.children { + holes[i] = ch.cnt + grands = append(grands, ch.children...) + } + out(bn.cnt, holes) + + // the grand children are new polygons. + for _, grand := range grands { + grand.generate(out) + } +} + +func (bn *bboxNode) generateRoot(out bboxOutFunc) { + for _, r := range bn.children { + r.generate(out) + } +} + +func buildMultipolygon( + cnts []contourmap.Contour, + reprojX, reprojY func(float64) float64, +) wkb.MultiPolygonGeom { + + var forest bboxNode + + for _, cnt := range cnts { + forest.insertRoot(contour(cnt)) + } + + //log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots)) + + var mp wkb.MultiPolygonGeom + + out := func(sh contour, hls []contour) { + + polygon := make(wkb.PolygonGeom, 1+len(hls)) + + // Handle shell + shell := make(wkb.LinearRingGeom, len(sh)) + for j, pt := range sh { + shell[j] = wkb.PointGeom{ + X: reprojX(pt.X), + Y: reprojY(pt.Y), + } + } + if shell.CCW() { + shell.Reverse() + } + polygon[0] = shell + + // handle holes + for i, hl := range hls { + hole := make(wkb.LinearRingGeom, len(hl)) + for j, pt := range hl { + hole[j] = wkb.PointGeom{ + X: reprojX(pt.X), + Y: reprojY(pt.Y), + } + } + if !hole.CCW() { + hole.Reverse() + } + polygon[1+i] = hole + } + + mp = append(mp, polygon) + } + + forest.generateRoot(out) + + return mp +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/strtree.go Tue Nov 05 14:30:22 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 mesh + +import ( + "bytes" + "compress/gzip" + "encoding/binary" + "io" + "log" + "math" + "sort" +) + +const STRTreeDefaultEntries = 8 + +type STRTree struct { + Entries int + tin *Tin + index []int32 + 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 + + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + + if top > 0 { // node + if s.bbox(top).Contains(x, y) { + n := s.index[top+1] + stack = append(stack, s.index[top+2:top+2+n]...) + } + } else { // leaf + top = -top - 1 + if !s.bbox(top).Contains(x, y) { + 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]], + } + if t.Contains(x, y) { + return t.Plane3D().Z(x, y), true + } + } + } + } + return 0, false +} + +func (s *STRTree) Build(t *Tin) { + + if s.Entries == 0 { + s.Entries = STRTreeDefaultEntries + } + + s.tin = t + + all := make([]int32, len(t.Triangles)) + + for i := range all { + all[i] = int32(i) + } + + s.index = append(s.index, 0) + + root := s.build(all) + + s.index[0] = root +} + +func (s *STRTree) BuildWithout(t *Tin, remove map[int32]struct{}) { + + if s.Entries == 0 { + s.Entries = STRTreeDefaultEntries + } + + s.tin = t + + all := make([]int32, 0, len(t.Triangles)-len(remove)) + + for i := 0; i < len(t.Triangles); i++ { + idx := int32(i) + if _, found := remove[idx]; !found { + all = append(all, idx) + } + } + + s.index = append(s.index, 0) + + root := s.build(all) + + s.index[0] = root +} + +func (s *STRTree) Clip(p *Polygon) map[int32]struct{} { + + removed := make(map[int32]struct{}) + + stack := []int32{s.index[0]} + + vertices := s.tin.Vertices + + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + + if top > 0 { // node + switch p.IntersectionBox2D(s.bbox(top)) { + case IntersectionInside: + // all triangles are inside polygon + case IntersectionOutSide: + // all triangles are outside polygon + s.allTriangles(top, removed) + default: + // mixed bag + n := s.index[top+1] + stack = append(stack, s.index[top+2:top+2+n]...) + } + } else { // leaf + top = -top - 1 + switch p.IntersectionBox2D(s.bbox(top)) { + case IntersectionInside: + // all triangles are inside polygon + case IntersectionOutSide: + // all triangles are outside polygon + n := s.index[top+1] + for _, idx := range s.index[top+2 : top+2+n] { + removed[idx] = struct{}{} + } + default: + 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]], + } + if p.IntersectionWithTriangle(&t) != IntersectionInside { + removed[idx] = struct{}{} + } + } + } + } + } + + 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 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + if top > 0 { // node + n := s.index[top+1] + stack = append(stack, s.index[top+2:top+2+n]...) + } else { // leaf + top = -top - 1 + n := s.index[top+1] + for _, idx := range s.index[top+2 : top+2+n] { + tris[idx] = struct{}{} + } + } + } +} + +func (s *STRTree) build(items []int32) int32 { + sort.Slice(items, func(i, j int) bool { + ti := s.tin.Triangles[items[i]] + tj := s.tin.Triangles[items[j]] + return s.tin.Vertices[ti[0]].X < s.tin.Vertices[tj[0]].X + }) + + P := int(math.Ceil(float64(len(items)) / float64(s.Entries))) + S := int(math.Ceil(math.Sqrt(float64(P)))) + + slices := s.strSplit(items, S) + + leaves := s.strJoin( + slices, S, + func(i, j int32) bool { + ti := s.tin.Triangles[i] + tj := s.tin.Triangles[j] + return s.tin.Vertices[ti[0]].Y < s.tin.Vertices[tj[0]].Y + }, + s.allocLeaf, + ) + + return s.buildNodes(leaves) +} + +func (s *STRTree) buildNodes(items []int32) int32 { + + if len(items) <= s.Entries { + return s.allocNode(items) + } + + sort.Slice(items, func(i, j int) bool { + return s.bbox(items[i]).X1 < s.bbox(items[j]).X1 + }) + + P := int(math.Ceil(float64(len(items)) / float64(s.Entries))) + S := int(math.Ceil(math.Sqrt(float64(P)))) + + slices := s.strSplit(items, S) + + nodes := s.strJoin( + slices, S, + func(i, j int32) bool { return s.bbox(i).Y1 < s.bbox(j).Y1 }, + s.allocNode, + ) + + return s.buildNodes(nodes) +} + +func (s *STRTree) bbox(idx int32) Box2D { + if idx < 0 { // Don't care if leaf or node. + idx = -idx - 1 + } + return s.bboxes[s.index[idx]] +} + +func (s *STRTree) strSplit(items []int32, S int) [][]int32 { + sm := S * s.Entries + slices := make([][]int32, S) + for i := range slices { + var n int + if l := len(items); l < sm { + n = l + } else { + n = sm + } + slices[i] = items[:n] + items = items[n:] + } + return slices +} + +func (s *STRTree) strJoin( + slices [][]int32, S int, + less func(int32, int32) bool, + alloc func([]int32) int32, +) []int32 { + nodes := make([]int32, 0, S*S) + + for _, slice := range slices { + sort.Slice(slice, func(i, j int) bool { + return less(slice[i], slice[j]) + }) + + for len(slice) > 0 { + var n int + if l := len(slice); l >= s.Entries { + n = s.Entries + } else { + n = l + } + nodes = append(nodes, alloc(slice[:n])) + slice = slice[n:] + } + } + return nodes +} + +func (s *STRTree) allocNode(items []int32) int32 { + pos := len(s.index) + s.index = append(s.index, 0, int32(len(items))) + s.index = append(s.index, items...) + if len(items) > 0 { + box := s.bbox(items[0]) + for i := 1; i < len(items); i++ { + box = box.Union(s.bbox(items[i])) + } + s.index[pos] = int32(s.allocBBox(box)) + } + return int32(pos) +} + +func (s *STRTree) allocBBox(box Box2D) int { + pos := len(s.bboxes) + s.bboxes = append(s.bboxes, box) + return pos +} + +func (s *STRTree) allocLeaf(items []int32) int32 { + pos := len(s.index) + s.index = append(s.index, 0, int32(len(items))) + s.index = append(s.index, items...) + if len(items) > 0 { + vertices := s.tin.Vertices + ti := s.tin.Triangles[items[0]] + t := Triangle{ + vertices[ti[0]], + vertices[ti[1]], + vertices[ti[2]], + } + box := t.BBox() + for i := 1; i < len(items); i++ { + it := items[i] + ti := s.tin.Triangles[it] + t := Triangle{ + vertices[ti[0]], + vertices[ti[1]], + vertices[ti[2]], + } + box = box.Union(t.BBox()) + } + s.index[pos] = int32(s.allocBBox(box)) + } + return int32(-(pos + 1)) +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/tin.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,271 @@ +// 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 mesh + +import ( + "bytes" + "encoding/binary" + "errors" + "fmt" + "io" + "log" + "math" + + "gemma.intevation.de/gemma/pkg/models" + "gemma.intevation.de/gemma/pkg/wkb" +) + +var ( + errNoByteSlice = errors.New("Not a byte slice") + errTooLessPoints = errors.New("Too less points") +) + +// Tin stores a mesh of triangles with common vertices. +type Tin struct { + // EPSG holds the projection. + EPSG uint32 + // Vertices are the shared vertices. + Vertices []Vertex + // Triangles are the triangles. + Triangles [][]int32 + + // Min is the lower left corner of the bbox. + Min Vertex + // Max is the upper right corner of the bbox. + Max Vertex +} + +func (t *Tin) Clip(polygon *Polygon) map[int32]struct{} { + var tree STRTree + tree.Build(t) + return tree.Clip(polygon) +} + +// FromWKB constructs the TIN from a WKB representation. +// Shared vertices are identified and referenced by the +// same index. +func (t *Tin) FromWKB(data []byte) error { + log.Printf("info: data length %d\n", len(data)) + + 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.TinZ: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var num uint32 + if err = binary.Read(r, order, &num); err != nil { + return err + } + + vertices := make([]Vertex, 0, 100000) + + var v Vertex + + v2i := make(map[Vertex]int32, 100000) + + var indexPool []int32 + + allocIndices := func() []int32 { + if len(indexPool) == 0 { + indexPool = make([]int32, 3*8*1024) + } + ids := indexPool[:3] + indexPool = indexPool[3:] + return ids + } + + var triangles [][]int32 + + min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} + max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} + + for i := uint32(0); i < num; i++ { + + endian, err = r.ReadByte() + 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) + } + + err = binary.Read(r, order, &geomType) + switch { + case err != nil: + return err + case geomType != wkb.TriangleZ: + return fmt.Errorf("unknown geometry type %d", geomType) + } + + var rings uint32 + if err = binary.Read(r, order, &rings); err != nil { + return err + } + triangle := allocIndices() + + for ring := uint32(0); ring < rings; ring++ { + var npoints uint32 + if err = binary.Read(r, order, &npoints); err != nil { + return err + } + + if npoints < 3 { + return errTooLessPoints + } + + for p := uint32(0); p < npoints; p++ { + var x, y, z uint64 + for _, addr := range []*uint64{&x, &y, &z} { + if err = binary.Read(r, order, addr); err != nil { + return err + } + } + if p >= 3 || ring >= 1 { + // Don't store the forth point. + continue + } + // Do this conversion later to spare reflect calls + // and allocs in binary.Read. + v.X = math.Float64frombits(x) + v.Y = math.Float64frombits(y) + v.Z = math.Float64frombits(z) + idx, found := v2i[v] + if !found { + idx = int32(len(vertices)) + v2i[v] = idx + vertices = append(vertices, v) + min.Minimize(v) + max.Maximize(v) + } + triangle[p] = idx + } + } + triangles = append(triangles, triangle) + } + + log.Printf("info: bbox: [[%f, %f], [%f, %f]]\n", + min.X, min.Y, max.X, max.Y) + + *t = Tin{ + EPSG: models.WGS84, + Vertices: vertices, + Triangles: triangles, + Min: min, + Max: max, + } + + return nil +} + +// Scan implements the sql.Scanner interface. +func (t *Tin) Scan(raw interface{}) error { + if raw == nil { + return nil + } + data, ok := raw.([]byte) + if !ok { + return errNoByteSlice + } + return t.FromWKB(data) +} + +func (t *Tin) serialize(w io.Writer) error { + + if err := binary.Write(w, binary.LittleEndian, t.EPSG); err != nil { + return err + } + + if err := t.Min.Write(w); err != nil { + return err + } + if err := t.Max.Write(w); err != nil { + return err + } + + if err := binary.Write( + w, binary.LittleEndian, uint32(len(t.Vertices))); err != nil { + 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 { + 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) + + if err := binary.Write( + w, binary.LittleEndian, uint32(len(t.Triangles))); err != nil { + return err + } + + var buf [binary.MaxVarintLen32]byte + var written int + var last int32 + for _, triangle := range t.Triangles { + for _, idx := range triangle { + value := idx - last + n := binary.PutVarint(buf[:], int64(value)) + 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 = idx + } + } + log.Printf("info: compressed tin indices in bytes: %d (%d)\n", + written, 3*4*len(t.Triangles)) + + return nil +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/triangulation.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,318 @@ +// 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 mesh + +import ( + "fmt" + "log" + "math" + + "gonum.org/v1/gonum/stat" +) + +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) EstimateTooLong() float64 { + + num := len(t.Triangles) / 3 + + lengths := make([]float64, 0, num) + + points := t.Points + +tris: + for i := 0; i < num; i++ { + idx := i * 3 + var max float64 + vs := t.Triangles[idx : idx+3] + for j, vj := range vs { + if t.Halfedges[idx+j] < 0 { + continue tris + } + k := (j + 1) % 3 + if l := points[vj].Distance2D(points[vs[k]]); l > max { + max = l + } + } + lengths = append(lengths, max) + } + + std := stat.StdDev(lengths, nil) + return 3.5 * std +} + +func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) { + + if tooLong <= 0 { + tooLong = t.EstimateTooLong() + } + + tooLong *= tooLong + + var candidates []int32 + + points := t.Points + + for i, num := 0, len(t.Triangles)/3; i < num; i++ { + idx := i * 3 + var max float64 + vs := t.Triangles[idx : idx+3] + for j, vj := range vs { + k := (j + 1) % 3 + if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max { + max = l + } + } + if max > tooLong { + candidates = append(candidates, int32(i)) + } + } + + removed := map[int32]struct{}{} + + isBorder := func(n int32) bool { + n *= 3 + for i := int32(0); i < 3; i++ { + e := n + i + o := t.Halfedges[e] + if o < 0 { + return true + } + if _, found := removed[o/3]; found { + return true + } + } + return false + } + + var newCandidates []int32 + + log.Printf("info: candidates: %d\n", len(candidates)) + for len(candidates) > 0 { + + oldRemoved := len(removed) + + for _, i := range candidates { + + if isBorder(i) { + removed[i] = struct{}{} + } else { + newCandidates = append(newCandidates, i) + } + } + + if oldRemoved == len(removed) { + break + } + + candidates = newCandidates + newCandidates = newCandidates[:0] + } + + log.Printf("info: candidates left: %d\n", len(candidates)) + log.Printf("info: triangles: %d\n", len(t.Triangles)/3) + log.Printf("info: triangles to remove: %d\n", len(removed)) + + type edge struct { + a, b int32 + prev, next *edge + } + + isClosed := func(e *edge) bool { + for curr := e.next; curr != nil; curr = curr.next { + if curr == e { + return true + } + } + return false + } + + open := map[int32]*edge{} + var rings []*edge + + for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ { + if _, found := removed[i]; found { + continue + } + n := i * 3 + for j := int32(0); j < 3; j++ { + e := n + j + f := t.Halfedges[e] + if f >= 0 { + if _, found := removed[f/3]; !found { + continue + } + } + a := t.Triangles[e] + b := t.Triangles[n+(j+1)%3] + + curr := &edge{a: a, b: b} + + if old := open[a]; old != nil { + delete(open, a) + if old.a == a { + old.prev = curr + curr.next = old + } else { + old.next = curr + curr.prev = old + } + + if isClosed(old) { + rings = append(rings, old) + } + } else { + open[a] = curr + } + + if old := open[b]; old != nil { + delete(open, b) + if old.b == b { + old.next = curr + curr.prev = old + } else { + old.prev = curr + curr.next = old + } + + if isClosed(old) { + rings = append(rings, old) + } + } else { + open[b] = curr + } + } + } + + if len(open) > 0 { + log.Printf("warn: open vertices left: %d\n", len(open)) + } + + if len(rings) == 0 { + log.Println("warn: no ring found") + return nil, removed + } + + curr := rings[0] + + polygon := LineStringZ{ + points[curr.a], + points[curr.b], + } + + for curr = curr.next; curr != rings[0]; curr = curr.next { + polygon = append(polygon, points[curr.b]) + } + + polygon = append(polygon, t.Points[rings[0].a]) + + log.Printf("length of boundary: %d\n", len(polygon)) + + return polygon, removed +} + +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/mesh/triangulator.go Tue Nov 05 14:30:22 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 mesh + +import ( + "errors" + "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 (tri *triangulator) Len() int { + return len(tri.points) +} + +func (tri *triangulator) Swap(i, j int) { + tri.ids[i], tri.ids[j] = tri.ids[j], tri.ids[i] +} + +func (tri *triangulator) Less(i, j int) bool { + d1 := tri.squaredDistances[tri.ids[i]] + d2 := tri.squaredDistances[tri.ids[j]] + if d1 != d2 { + return d1 < d2 + } + p1 := tri.points[tri.ids[i]] + p2 := tri.points[tri.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 errors.New("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 (tri *triangulator) hashKey(point Vertex) int { + d := point.Sub2D(tri.center) + return int(pseudoAngle(d.X, d.Y) * float64(len(tri.hash))) +} + +func (tri *triangulator) hashEdge(e *node) { + tri.hash[tri.hashKey(tri.points[e.i])] = e +} + +func (tri *triangulator) addTriangle(i0, i1, i2, a, b, c int32) int32 { + i := int32(tri.trianglesLen) + tri.triangles[i] = i0 + tri.triangles[i+1] = i1 + tri.triangles[i+2] = i2 + tri.link(i, a) + tri.link(i+1, b) + tri.link(i+2, c) + tri.trianglesLen += 3 + return i +} + +func (tri *triangulator) link(a, b int32) { + tri.halfedges[a] = b + if b >= 0 { + tri.halfedges[b] = a + } +} + +func (tri *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 := tri.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 := tri.triangles[ar] + pr := tri.triangles[a] + pl := tri.triangles[al] + p1 := tri.triangles[bl] + + illegal := inCircle(tri.points[p0], tri.points[pr], tri.points[pl], tri.points[p1]) + + if illegal { + tri.triangles[a] = p1 + tri.triangles[b] = p0 + + // edge swapped on the other side of the hull (rare) + // fix the halfedge reference + if tri.halfedges[bl] == -1 { + e := tri.hull + for { + if e.t == bl { + e.t = a + break + } + e = e.next + if e == tri.hull { + break + } + } + } + + tri.link(a, tri.halfedges[bl]) + tri.link(b, tri.halfedges[ar]) + tri.link(ar, bl) + + br := b0 + (b+1)%3 + + tri.legalize(a) + return tri.legalize(br) + } + + return ar +} + +func (tri *triangulator) convexHull() []Vertex { + var result []Vertex + e := tri.hull + for e != nil { + result = append(result, tri.points[e.i]) + e = e.prev + if e == tri.hull { + break + } + } + return result +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/util.go Tue Nov 05 14:30:22 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 mesh + +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)) +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/vertex.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,1229 @@ +// 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 mesh + +import ( + "bytes" + "encoding/binary" + "fmt" + "io" + "log" + "math" + "sort" + + "gemma.intevation.de/gemma/pkg/wkb" +) + +type ( + Point struct { + X float64 + Y float64 + } + + // Vertex is a 3D vertex. + Vertex struct { + X float64 + Y float64 + Z float64 + } + + // Triangle is a triangle consisting of three vertices. + Triangle [3]Vertex + + // Line is a line defined by first vertex on that line + // and the second being the direction. + Line [2]Vertex + + // Box is a 3D box. + Box [2]Vertex + + // MultiPointZ is a set of vertices. + MultiPointZ []Vertex + + // LineStringZ is a line string formed of vertices. + LineStringZ []Vertex + + // MultiLineStringZ is a set of line strings. + MultiLineStringZ []LineStringZ + + // Box2D is 2D area from (X1, Y1) to (X2, Y2). + Box2D struct { + X1 float64 + Y1 float64 + X2 float64 + Y2 float64 + } + + // Plane2D is a 2D plane (a line in 2D space). + Plane2D struct { + A float64 + 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 (t *Triangle) BBox() Box2D { + minX := math.Min(math.Min(t[0].X, t[1].X), t[2].X) + maxX := math.Max(math.Max(t[0].X, t[1].X), t[2].X) + minY := math.Min(math.Min(t[0].Y, t[1].Y), t[2].Y) + maxY := math.Max(math.Max(t[0].Y, t[1].Y), t[2].Y) + return Box2D{ + X1: minX, Y1: minY, + X2: maxX, Y2: maxY, + } +} + +func (a Box2D) Inside(b Box2D) bool { + return a.X1 >= b.X1 && a.X2 <= b.X2 && + a.Y1 >= b.Y1 && a.Y2 <= b.Y2 +} + +func (a Box2D) Size() (float64, float64) { + return a.X2 - a.X1, a.Y2 - a.Y1 +} + +func (a Box2D) Empty() bool { + const eps = 0.0000001 + return math.Abs(a.X2-a.X1) < eps && + math.Abs(a.Y2-a.Y1) < eps +} + +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 (p Plane3D) Eval(v Vertex) float64 { + return p.A*v.X + p.B*v.Y + p.C*v.Z + p.D +} + +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) +} + +func (v Vertex) Distance(w Vertex) float64 { + v = v.Sub(w) + 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) { + if w.X < v.X { + v.X = w.X + } + if w.Y < v.Y { + v.Y = w.Y + } + if w.Z < v.Z { + v.Z = w.Z + } +} + +// Maximize adjust this vertex v to hold the maximum +// values component-wise of v and w. +func (v *Vertex) Maximize(w Vertex) { + if w.X > v.X { + v.X = w.X + } + if w.Y > v.Y { + v.Y = w.Y + } + if w.Z > v.Z { + v.Z = w.Z + } +} + +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{ + v.X - w.X, + v.Y - w.Y, + v.Z - w.Z, + } +} + +// Add returns (v + w) component-wise. +func (v Vertex) Add(w Vertex) Vertex { + return Vertex{ + v.X + w.X, + v.Y + w.Y, + v.Z + w.Z, + } +} + +// Scale returns s*v component-wise. +func (v Vertex) Scale(s float64) Vertex { + return Vertex{ + s * v.X, + s * v.Y, + s * v.Z, + } +} + +// Interpolate returns a function that return s*b[1] + b[0] +// component-wise. +func (b Box) Interpolate() func(Vertex) Vertex { + v1, v2 := b[0], b[1] + v2 = v2.Sub(v1) + return func(s Vertex) Vertex { + return Vertex{ + v2.X*s.X + v1.X, + v2.Y*s.Y + v1.Y, + v2.Z*s.Z + v1.Z, + } + } +} + +func (b Box) HasX() bool { return math.Abs(b[0].X-b[1].X) > epsPlane } +func (b Box) HasY() bool { return math.Abs(b[0].Y-b[1].Y) > epsPlane } +func (b Box) HasZ() bool { return math.Abs(b[0].Z-b[1].Z) > epsPlane } + +// Less returns if one of v component is less than the +// corresponing component in w. +func (v Vertex) Less(w Vertex) bool { + return v.X < w.X || v.Y < w.Y || v.Z < w.Z +} + +// NewLine return a line of point/direction. +func NewLine(p1, p2 Vertex) Line { + return Line{ + p2.Sub(p1), + p1, + } +} + +// Eval returns the vertex for t*l[0] + l[1]. +func (l Line) Eval(t float64) Vertex { + return l[0].Scale(t).Add(l[1]) +} + +// IntersectHorizontal returns the intersection point +// for a given z value. +func (l Line) IntersectHorizontal(h float64) Vertex { + t := (h - l[1].Z) / l[0].Z + return l.Eval(t) +} + +func side(z, h float64) int { + switch { + case z < h: + return -1 + case z > h: + return +1 + } + 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].Sub2D(t[0]) + v1 := t[1].Sub2D(t[0]) + v2 := Vertex{X: x, Y: y}.Sub2D(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), +// one vertex (only touching an vertex) or +// two vertices (real intersection). +func (t *Triangle) IntersectHorizontal(h float64) LineStringZ { + sides := [3]int{ + side(t[0].Z, h), + side(t[1].Z, h), + side(t[2].Z, h), + } + + var points LineStringZ + + for i := 0; i < 3; i++ { + j := (i + 1) % 3 + si := sides[i] + sj := sides[j] + + switch { + case si == 0 && sj == 0: + // both on plane + points = append(points, t[i], t[j]) + case si == 0 && sj != 0: + // first on plane + points = append(points, t[i]) + case si != 0 && sj == 0: + // second on plane + points = append(points, t[j]) + case si == sj: + // both on same side + default: + // real intersection + v := NewLine(t[i], t[j]).IntersectHorizontal(h) + points = append(points, v) + } + } + + return points +} + +func linearScale(x1, y1, x2, y2 float64) func(Vertex) float64 { + dx := x2 - x1 + dy := y2 - y1 + + switch { + case dx != 0: + return func(v Vertex) float64 { + return (v.X - x1) / dx + } + case dy != 0: + return func(v Vertex) float64 { + return (v.Y - y1) / dy + } + default: + return func(Vertex) float64 { + return 0 + } + } +} + +func (ls LineStringZ) BBox() Box2D { + + min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} + max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} + + for _, v := range ls { + min.Minimize(v) + max.Maximize(v) + } + + return Box2D{ + X1: min.X, + Y1: min.Y, + X2: max.X, + Y2: max.Y, + } +} + +func (ls LineStringZ) Area() float64 { + return polygonArea(ls) +} + +func (ls LineStringZ) Reverse() { + for i, j := 0, len(ls)-1; i < j; i, j = i+1, j-1 { + ls[i], ls[j] = ls[j], ls[i] + } +} + +func (ls LineStringZ) order(position func(Vertex) float64) { + type posVertex struct { + pos float64 + v Vertex + } + positions := make([]posVertex, len(ls)) + for i, v := range ls { + positions[i] = posVertex{position(v), v} + } + sort.Slice(positions, func(i, j int) bool { + return positions[i].pos < positions[j].pos + }) + for i := range positions { + ls[i] = positions[i].v + } +} + +// EpsEquals returns true if v and w are equal component-wise +// with the values within a epsilon range. +func (v Vertex) EpsEquals(w Vertex) bool { + const eps = 1e-5 + return math.Abs(v.X-w.X) < eps && + math.Abs(v.Y-w.Y) < eps && math.Abs(v.Z-w.Z) < eps +} + +// EpsEquals returns true if v and w are equal component-wise +// in the X/Y plane with the values within a epsilon range. +func (v Vertex) EpsEquals2D(w Vertex) bool { + const eps = 1e-5 + return math.Abs(v.X-w.X) < eps && + math.Abs(v.Y-w.Y) < eps +} + +// JoinOnLine joins the the elements of a given multi line string +// under the assumption that the segments are all on the line segment +// from (x1, y1) to (x2, y2). +func (mls MultiLineStringZ) JoinOnLine(x1, y1, x2, y2 float64) MultiLineStringZ { + + position := linearScale(x1, y1, x2, y2) + + type posLineString struct { + pos float64 + line LineStringZ + } + + positions := make([]posLineString, 0, len(mls)) + + for _, ls := range mls { + if len(ls) == 0 { + continue + } + // order the atoms first + ls.order(position) + positions = append(positions, posLineString{position(ls[0]), ls}) + } + + sort.Slice(positions, func(i, j int) bool { + return positions[i].pos < positions[j].pos + }) + + out := make(MultiLineStringZ, 0, len(positions)) + + var ignored int + + for i := range positions { + curr := positions[i].line + if l := len(out); l > 0 { + last := out[l-1] + + if last[len(last)-1].EpsEquals(curr[0]) { + out[l-1] = append(last[:len(last)-1], curr...) + continue + } + if position(last[len(last)-1]) > position(curr[0]) { + ignored++ + continue + } + } + out = append(out, curr) + } + + log.Printf("info: ignored parts: %d\n", ignored) + + return out +} + +// Write writes a Vertex as three 64 bit values in little endian order +// to the given writer. +func (v *Vertex) Write(w io.Writer) error { + if err := binary.Write( + w, binary.LittleEndian, math.Float64bits(v.X)); err != nil { + return err + } + if err := binary.Write( + w, binary.LittleEndian, math.Float64bits(v.Y)); err != nil { + return err + } + return binary.Write( + w, binary.LittleEndian, math.Float64bits(v.Z)) +} + +// Read fills this vertex with three 64 bit values stored as +// little endian from the given reader. +func (v *Vertex) Read(r io.Reader) error { + var buf [8]byte + b := buf[:] + if _, err := io.ReadFull(r, b); err != nil { + return nil + } + v.X = math.Float64frombits(binary.LittleEndian.Uint64(b)) + if _, err := io.ReadFull(r, b); err != nil { + return nil + } + v.Y = math.Float64frombits(binary.LittleEndian.Uint64(b)) + if _, err := io.ReadFull(r, b); err != nil { + return nil + } + v.Z = math.Float64frombits(binary.LittleEndian.Uint64(b)) + return nil +} + +// AsWKB returns the WKB representation of the given multi line string. +func (mls MultiLineStringZ) AsWKB() []byte { + + // pre-calculate size to avoid reallocations. + size := 1 + 4 + 4 + for _, ml := range mls { + size += 1 + 4 + 4 + len(ml)*3*8 + } + + buf := bytes.NewBuffer(make([]byte, 0, size)) + + 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, 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)) + binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) + binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z)) + } + } + + return buf.Bytes() +} + +// AsWKB2D returns the WKB representation of the given multi line string +// leaving the z component out. +func (mls MultiLineStringZ) AsWKB2D() []byte { + + // pre-calculate size to avoid reallocations. + size := 1 + 4 + 4 + for _, ml := range mls { + size += 1 + 4 + 4 + len(ml)*2*8 + } + + buf := bytes.NewBuffer(make([]byte, 0, size)) + + 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, 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)) + binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) + } + } + + return buf.Bytes() +} + +func (ls LineStringZ) CCW() bool { + var sum float64 + for i, v1 := range ls { + v2 := ls[(i+1)%len(ls)] + sum += (v2.X - v1.X) * (v2.Y + v1.Y) + } + return sum > 0 +} + +// Join joins two lines leaving the first of the second out. +func (ls LineStringZ) Join(other LineStringZ) LineStringZ { + nline := make(LineStringZ, len(ls)+len(other)-1) + copy(nline, ls) + copy(nline[len(ls):], other[1:]) + return nline +} + +// Merge merges line segments of a given multi line string +// by finding common start and end vertices. +func (mls MultiLineStringZ) Merge() MultiLineStringZ { + + var out MultiLineStringZ + + min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} + + for _, line := range mls { + for _, v := range line { + min.Minimize(v) + } + } + + type point struct { + x int64 + y int64 + } + + const precision = 1e7 + + quant := func(v Vertex) point { + return point{ + x: int64(math.Round((v.X - min.X) * precision)), + y: int64(math.Round((v.Y - min.Y) * precision)), + } + } + + heads := make(map[point]*[]LineStringZ) + + for _, line := range mls { + if len(line) < 2 { + out = append(out, line) + continue + } + head := quant(line[0]) + tail := quant(line[len(line)-1]) + if head == tail { // its already a ring + out = append(out, line) + continue + } + + if hs := heads[tail]; hs != nil { + l := len(*hs) + last := (*hs)[l-1] + if l == 1 { + delete(heads, tail) + } else { + (*hs)[l-1] = nil + *hs = (*hs)[:l-1] + } + line = line.Join(last) + + if head == quant(line[len(line)-1]) { // its a ring + out = append(out, line) + continue + } + } + + if hs := heads[head]; hs != nil { + *hs = append(*hs, line) + } else { + heads[head] = &[]LineStringZ{line} + } + } + +again: + for head, lines := range heads { + for i, line := range *lines { + tail := quant(line[len(line)-1]) + for hs := heads[tail]; hs != nil && len(*hs) > 0; hs = heads[tail] { + l := len(*hs) + last := (*hs)[l-1] + (*hs)[l-1] = nil + *hs = (*hs)[:l-1] + line = line.Join(last) + + if tail = quant(line[len(line)-1]); head == tail { // its a ring + out = append(out, line) + // remove from current lines + copy((*lines)[i:], (*lines)[i+1:]) + (*lines)[len(*lines)-1] = nil + *lines = (*lines)[:len(*lines)-1] + goto again + } + // overwrite in current lines + (*lines)[i] = line + } + } + } + + // rings := len(out) + + // The rest are open line strings. + for _, lines := range heads { + for _, line := range *lines { + out = append(out, line) + } + } + + // log.Printf("segments before/after merge: %d/%d (%d rings)\n", + // len(mls), len(out), rings) + + 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 { + return a.X1 + } + return a.X2 +} + +// Yi returns the i-th y component. +func (a Box2D) Yi(i int) float64 { + if i == 0 { + return a.Y1 + } + return a.Y2 +} + +func (a Box2D) Union(b Box2D) Box2D { + return Box2D{ + X1: math.Min(a.X1, b.X1), + Y1: math.Min(a.Y1, b.Y1), + X2: math.Max(a.X2, b.X2), + Y2: math.Max(a.Y2, b.Y2), + } +} + +func (a Box2D) Area() float64 { + return (a.X2 - a.X1) * (a.Y2 - a.Y1) +} + +// NewPlane2D creates a new Plane2D from two given points. +func NewPlane2D(x1, y1, x2, y2 float64) Plane2D { + b := x2 - x1 + a := -(y2 - y1) + + l := math.Sqrt(a*a + b*b) + a /= l + b /= l + + // a*x1 + b*y1 + c = 0 + c := -(a*x1 + b*y1) + return Plane2D{a, b, c} +} + +// Eval determines the distance of a given point +// from the plane. The sign of the result indicates +// the sideness. +func (p Plane2D) Eval(x, y float64) float64 { + return p.A*x + p.B*y + p.C +} + +const epsPlane = 1e-5 + +func sides(s int, x float64) int { + if math.Signbit(x) { + return s | 2 + } + return s | 1 +} + +// IntersectsPlane checks if a Box2D intersects with +// a given Plane2D. +func (a Box2D) IntersectsPlane(p Plane2D) bool { + var s int + for i := 0; i < 2; i++ { + x := a.Xi(i) + for j := 0; j < 2; j++ { + y := a.Yi(j) + v := p.Eval(x, y) + if math.Abs(v) < epsPlane { + //log.Printf("on line") + return true + } + if s = sides(s, v); s == 3 { + //log.Printf("... on both sides (djt)") + return true + } + } + } + //log.Printf("side: %d\n", s) + return false +} + +// Cross calculates the cross product of two vertices. +func (v Vertex) Cross(w Vertex) Vertex { + return Vertex{ + v.Y*w.Z - v.Z*w.Y, + v.Z*w.X - v.X*w.Z, + v.X*w.Y - v.Y*w.X, + } +} + +// Intersection calcultes the 2D intersection point of +// two Plane2Ds. If they do not intersect the returned +// bool flags is set to false. +func (p Plane2D) Intersection(o Plane2D) (float64, float64, bool) { + + u1 := Vertex{p.A, p.B, p.C} + u2 := Vertex{o.A, o.B, o.C} + + plane := u1.Cross(u2) + + if plane.Z == 0 { + return 0, 0, false + } + + return plane.X / plane.Z, plane.Y / plane.Z, true +} + +// VerticalLine is a 2D line segment. +type VerticalLine struct { + x1 float64 + y1 float64 + x2 float64 + y2 float64 + + line Plane2D +} + +// NewVerticalLine creates a new 2D line segment. +func NewVerticalLine(x1, y1, x2, y2 float64) *VerticalLine { + return &VerticalLine{ + x1: x1, + y1: y1, + x2: x2, + y2: y2, + line: NewPlane2D(x1, y1, x2, y2), + } +} + +func onPlane(x float64) bool { return math.Abs(x) < epsPlane } + +func relative(v1, v2 Vertex) func(x, y float64) float64 { + ls := linearScale(v1.X, v1.Y, v2.X, v2.Y) + return func(x, y float64) float64 { + return ls(Vertex{x, y, 0}) + } +} + +func interpolate(a, b float64) func(float64) float64 { + return func(x float64) float64 { + return (b-a)*x + a + } +} + +func linear(v1, v2 Vertex) func(t float64) Vertex { + return func(t float64) Vertex { + return Vertex{ + (v2.X-v1.X)*t + v1.X, + (v2.Y-v1.Y)*t + v1.Y, + (v2.Z-v1.Z)*t + v1.Z, + } + } +} + +func inRange(a float64) bool { return 0 <= a && a <= 1 } + +// Intersection intersects a line segment with a triangle. +func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ { + + var out LineStringZ + + //defer func() { log.Printf("length out: %d\n", len(out)) }() + +edges: + for i := 0; i < 3 && len(out) < 2; i++ { + j := (i + 1) % 3 + edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y) + + s1 := edge.Eval(vl.x1, vl.y1) + s2 := edge.Eval(vl.x2, vl.y2) + + o1 := onPlane(s1) + o2 := onPlane(s2) + + // log.Printf("s1, s2: %t %t (%f %f)\n", o1, o2, s1, s2) + + switch { + case o1 && o2: + pos := relative(t[i], t[j]) + t1 := pos(vl.x1, vl.y1) + t2 := pos(vl.x2, vl.y2) + + r1 := inRange(t1) + r2 := inRange(t2) + + switch { + case r1 && r2: + lin := linear(t[i], t[j]) + out = append(out, lin(t1), lin(t2)) + return out + + case !r1 && !r2: // the hole edge + out = append(out, t[i], t[j]) + return out + case !r1: + if t1 < 0 { + // below first -> clip by first + lin := linear(t[i], t[j]) + out = append(out, t[i], lin(t2)) + } else { + // above second -> clip by second + lin := linear(t[i], t[j]) + out = append(out, lin(t2), t[j]) + } + return out + case !r2: + if t2 < 0 { + // below first -> clip by first + lin := linear(t[i], t[j]) + out = append(out, t[i], lin(t1)) + } else { + // above second -> clip by second + lin := linear(t[i], t[j]) + out = append(out, lin(t1), t[j]) + } + return out + } + + case o1: + t1 := relative(t[i], t[j])(vl.x1, vl.y1) + if !inRange(t1) { + continue edges + } + out = append(out, linear(t[i], t[j])(t1)) + + case o2: + t2 := relative(t[i], t[j])(vl.x2, vl.y2) + if !inRange(t2) { + continue edges + } + out = append(out, linear(t[i], t[j])(t2)) + + default: + x, y, intersects := vl.line.Intersection(edge) + if !intersects { + continue edges + } + + // log.Println("Intersection -----------------------------") + t1 := relative(t[i], t[j])(x, y) + // log.Printf("relative pos: %f\n", t1) + if !inRange(t1) { + continue edges + } + + // log.Println("valid ***************") + + z := interpolate(t[j].Z, t[i].Z)(t1) + n := Vertex{x, y, z} + + if math.Signbit(s1) != math.Signbit(s2) { + // log.Println("\ton different sides") + // different sides -> vertex on edge. + out = append(out, n) + } else { // same side -> inside. Find peer + if len(out) > 0 { // we already have our peer + out = append(out, n) + return out + } + + for k := 0; k < 3; k++ { + if k == i { + continue + } + l := (k + 1) % 3 + other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y) + xo, yo, intersects := vl.line.Intersection(other) + if !intersects { + continue + } + t2 := relative(t[k], t[l])(xo, yo) + if !inRange(t2) { + continue + } + zo := interpolate(t[k].Z, t[l].Z)(t2) + + m := Vertex{xo, yo, zo} + + var xn, yn, xf, yf float64 + + // Which is nearer to current edge? + if math.Abs(s1) < math.Abs(s2) { + xn, yn = vl.x1, vl.y1 + xf, yf = vl.x2, vl.y2 + } else { + xn, yn = vl.x2, vl.y2 + xf, yf = vl.x1, vl.y1 + } + + if onPlane(other.Eval(xn, yn)) { + // triangle intersect in only point + // treat as no intersection + break edges + } + + pos := relative(n, m) + + tzn := pos(xn, yn) + tzf := pos(xf, yf) + + if !inRange(tzn) { + // if near is out of range far is, too. + return out + } + + lin := interpolate(n.Z, m.Z) + + if inRange(tzf) { + m.X = xf + m.Y = yf + m.Z = lin(tzf) + } // else m is clipping + + n.X = xn + n.Y = yn + n.Z = lin(tzn) + + out = append(out, n, m) + return out + } + } + } + } + + // supress single point touches. + if len(out) == 1 { + out = out[:0] + } + + return out +} + +// AsWKB returns a WKB representation of the given point cloud. +func (mpz MultiPointZ) AsWKB() []byte { + size := 1 + 4 + 4 + len(mpz)*(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(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 { + 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)) + } + + return buf.Bytes() +} + +func (mpz *MultiPointZ) 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.MultiPointZ: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var numPoints uint32 + if err := binary.Read(r, order, &numPoints); err != nil { + return err + } + + points := make(MultiPointZ, numPoints) + + for i := range points { + endian, err = r.ReadByte() + + 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) + } + + err = binary.Read(r, order, &geomType) + + switch { + case err != nil: + return err + case geomType != wkb.PointZ: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var x, y, z uint64 + if err = binary.Read(r, order, &x); err != nil { + return err + } + if err = binary.Read(r, order, &y); err != nil { + return err + } + if err = binary.Read(r, order, &z); err != nil { + return err + } + points[i] = Vertex{ + X: math.Float64frombits(x), + Y: math.Float64frombits(y), + Z: math.Float64frombits(z), + } + } + + *mpz = points + + return nil +}
--- a/pkg/octree/areas.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,88 +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 - -import ( - "runtime" - "sync" - - "gemma.intevation.de/gemma/pkg/common" -) - -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 -}
--- a/pkg/octree/cache.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,199 +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 - -import ( - "context" - "database/sql" - "sync" - "time" -) - -type ( - cacheKey struct { - date time.Time - bottleneck string - } - - cacheEntry struct { - checksum string - tree *STRTree - access time.Time - } - - // Cache holds Octrees for a defined amount of time in memory - // before they are released. - Cache struct { - sync.Mutex - entries map[cacheKey]*cacheEntry - } -) - -const ( - cleanupCacheSleep = 6 * time.Minute - maxCacheAge = 5 * time.Minute - maxCacheEntries = 4 -) - -const ( - directMeshSQL = ` -SELECT mesh_index FROM waterway.sounding_results -WHERE id = $1 -` - fetchMeshSQL = ` -SELECT mesh_checksum, mesh_index -FROM waterway.sounding_results -WHERE bottleneck_id = $1 AND date_info = $2::date - AND mesh_checksum IS NOT NULL AND mesh_index IS NOT NULL -` - checkMeshSQL = ` -SELECT CASE - WHEN mesh_checksum = $3 THEN NULL - ELSE mesh_index - END -FROM waterway.sounding_results -WHERE bottleneck_id = $1 AND date_info = $2::date - AND mesh_checksum IS NOT NULL AND mesh_index IS NOT NULL -` -) - -var cache = Cache{ - entries: map[cacheKey]*cacheEntry{}, -} - -func init() { - go cache.background() -} - -func (c *Cache) background() { - for { - time.Sleep(cleanupCacheSleep) - c.cleanup() - } -} - -func (c *Cache) cleanup() { - c.Lock() - defer c.Unlock() - good := time.Now().Add(-maxCacheAge) - for k, v := range c.entries { - if v.access.Before(good) { - delete(c.entries, k) - } - } -} - -// FromCache fetches an Octree from the global Octree cache. -func FromCache( - ctx context.Context, - conn *sql.Conn, - bottleneck string, date time.Time, -) (*STRTree, error) { - return cache.get(ctx, conn, bottleneck, date) -} - -// FetchOctreeDirectly loads an octree directly from the database. -func FetchMeshDirectly( - ctx context.Context, - tx *sql.Tx, - id int64, -) (*STRTree, error) { - var data []byte - err := tx.QueryRowContext(ctx, directMeshSQL, id).Scan(&data) - if err != nil { - return nil, err - } - 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, -) (*STRTree, error) { - c.Lock() - defer c.Unlock() - - key := cacheKey{date, bottleneck} - entry := c.entries[key] - - var data []byte - var checksum string - - if entry == nil { - // fetch from database - err := conn.QueryRowContext( - ctx, fetchMeshSQL, bottleneck, date).Scan(&checksum, &data) - switch { - case err == sql.ErrNoRows: - return nil, nil - case err != nil: - return nil, err - } - } else { - // check if we are not outdated. - err := conn.QueryRowContext( - ctx, checkMeshSQL, bottleneck, date, entry.checksum).Scan(&data) - switch { - case err == sql.ErrNoRows: - return nil, nil - case err != nil: - return nil, err - } - if data == nil { // we are still current - entry.access = time.Now() - return entry.tree, nil - } - } - - tree := new(STRTree) - - if err := tree.FromBytes(data); err != nil { - return nil, err - } - - now := time.Now() - - if entry != nil { - entry.tree = tree - entry.access = now - return tree, nil - } - - for len(c.entries) >= maxCacheEntries { - // Evict the entry that is accessed the longest time ago. - var oldestKey cacheKey - oldest := now - - for k, v := range c.entries { - if v.access.Before(oldest) { - oldest = v.access - oldestKey = k - } - } - delete(c.entries, oldestKey) - } - - c.entries[key] = &cacheEntry{ - checksum: checksum, - tree: tree, - access: now, - } - - return tree, nil -}
--- a/pkg/octree/classbreaks.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,150 +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) 2019 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 ( - "context" - "database/sql" - "errors" - "math" - "sort" - "strconv" - "strings" -) - -const ( - selectClassBreaksSQL = ` -SELECT config_val FROM sys_admin.system_config -WHERE config_key = $1` -) - -func SampleDiffHeights(min, max, step float64) []float64 { - var heights []float64 - switch { - case min >= 0: // All values positive. - for v := 0.0; v <= max; v += step { - if v >= min { - heights = append(heights, v) - } - } - case max <= 0: // All values negative. - for v := 0.0; v >= min; v -= step { - if v <= max { - heights = append(heights, v) - } - } - default: // Positive and negative. - for v := step; v <= max; v += step { - 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 >= min; v -= step { - heights = append(heights, v) - } - } - return heights -} - -func ParseClassBreaks(config string) ([]float64, error) { - - parts := strings.Split(config, ",") - classes := make([]float64, 0, len(parts)) - for _, part := range parts { - if idx := strings.IndexRune(part, ':'); idx >= 0 { - part = part[:idx] - } - if part = strings.TrimSpace(part); part == "" { - continue - } - v, err := strconv.ParseFloat(part, 64) - if err != nil { - return nil, err - } - classes = append(classes, v) - } - - sort.Float64s(classes) - return classes, nil -} - -func LoadClassBreaks(ctx context.Context, tx *sql.Tx, key string) ([]float64, error) { - - var config sql.NullString - - err := tx.QueryRowContext(ctx, selectClassBreaksSQL, key).Scan(&config) - - switch { - case err == sql.ErrNoRows: - return nil, nil - case err != nil: - return nil, err - case !config.Valid: - return nil, errors.New("Invalid config string") - } - - return ParseClassBreaks(config.String) -} - -func round(v float64) float64 { - return math.Round(v*10000) / 10000 -} - -func ExtrapolateClassBreaks(cbs []float64, min, max float64) []float64 { - if min > max { - min, max = max, min - } - - n := make([]float64, len(cbs)) - copy(n, cbs) - sort.Float64s(n) - - for len(n) > 0 && n[0] < min { - n = n[1:] - } - - if len(n) == 0 { - return n - } - - for len(n) > 0 && n[len(n)-1] > max { - n = n[:len(n)-1] - } - - if len(n) == 0 { - return n - } - - for min < n[0] { - diff := n[1] - n[0] - if diff == 0 { - break - } - m := make([]float64, len(n)+1) - m[0] = round(n[0] - diff) - copy(m[1:], n) - n = m - } - - for max > n[len(n)-1] { - diff := n[len(n)-1] - n[len(n)-2] - if diff == 0 { - break - } - n = append(n, round(n[len(n)-1]+diff)) - } - - return n -}
--- a/pkg/octree/hull.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -// 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 Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,166 +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 - -import ( - "bufio" - "bytes" - "compress/gzip" - "encoding/binary" - "log" -) - -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 - } - - 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) - } - } - - for i := range bboxes { - read(&bboxes[i].X1) - read(&bboxes[i].Y1) - read(&bboxes[i].X2) - read(&bboxes[i].Y2) - } - - return 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 - } - - 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", - 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 err - } - - log.Printf("info: vertices: %d\n", numVertices) - - vertices := make([]Vertex, numVertices) - t.Vertices = vertices - - for i := range vertices { - if err := vertices[i].Read(r); err != nil { - return err - } - } - - var numTriangles uint32 - if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil { - return err - } - - log.Printf("info: triangles: %d\n", numTriangles) - - indices := make([]int32, 3*numTriangles) - triangles := make([][]int32, numTriangles) - t.Triangles = triangles - - var last int32 - - for i := range triangles { - tri := indices[:3] - indices = indices[3:] - triangles[i] = tri - for j := range tri { - v, err := binary.ReadVarint(r) - if err != nil { - return err - } - value := int32(v) + last - tri[j] = value - last = value - } - } - - return nil -}
--- a/pkg/octree/node.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,49 +0,0 @@ -// 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 -}
--- a/pkg/octree/plane2d_test.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +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 - -import ( - "math" - "testing" -) - -func TestIntersection(t *testing.T) { - - table := []struct { - a [4]float64 - b [4]float64 - intersects bool - x, y float64 - }{ - {[4]float64{-1, -1, 1, 1}, [4]float64{-1, 1, 1, -1}, true, 0, 0}, - {[4]float64{0, 0, 1, 1}, [4]float64{0, 1, 1, 0}, true, 0.5, 0.5}, - {[4]float64{0, 0, 1, 0}, [4]float64{0, 1, 1, 1}, false, 0, 0}, - } - - for _, e := range table { - p1 := NewPlane2D(e.a[0], e.a[1], e.a[2], e.a[3]) - p2 := NewPlane2D(e.b[0], e.b[1], e.b[2], e.b[3]) - x, y, intersects := p1.Intersection(p2) - if intersects != e.intersects { - t.Fatalf("Have %t want %t\n", intersects, e.intersects) - } - if e.intersects { - if math.Abs(e.x-x) > epsPlane || math.Abs(e.y-y) > epsPlane { - t.Fatalf("Have (%f, %f)t want (%f, %f)\n", x, y, e.x, e.y) - } - } - } -}
--- a/pkg/octree/polygon.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,505 +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 - -import ( - "bytes" - "encoding/binary" - "fmt" - "log" - "math" - - "github.com/tidwall/rtree" - - "gemma.intevation.de/gemma/pkg/wkb" -) - -type ( - ring []float64 - - Polygon struct { - 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 (ls lineSegment) intersects(a Box2D) 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 (ls lineSegment) intersectsLineSegment(o lineSegment) bool { - - p0 := ls[:2] - p1 := ls[2:4] - p2 := o[:2] - p3 := o[2:4] - - s10x := p1[0] - p0[0] - s10y := p1[1] - p0[1] - s32x := p3[0] - p2[0] - s32y := p3[1] - p2[1] - - den := s10x*s32y - s32x*s10y - - if den == 0 { - return false - } - - denPos := den > 0 - - s02x := p0[0] - p2[0] - s02y := p0[1] - p2[1] - - sNum := s10x*s02y - s10y*s02x - if sNum < 0 == denPos { - return false - } - - tNum := s32x*s02y - s32y*s02x - if tNum < 0 == denPos { - return false - } - - if sNum > den == denPos || tNum > den == denPos { - return false - } - - // t := tNum / den - // intersection at( p0[0] + (t * s10x), p0[1] + (t * s10y) ) - 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 { - if item.(lineSegment).intersects(box) { - intersects = true - return false - } - return true - }) - if intersects { - return IntersectionOverlaps - } - } - - // No intersection -> check inside or outside - // if an abritrary point is inside or not. - - // Check holes first: inside a hole means outside. - if len(p.rings) > 1 { - for _, hole := range p.rings[1:] { - if contains(hole, box.X1, box.Y1) { - return IntersectionOutSide - } - } - } - - // Check shell - if contains(p.rings[0], box.X1, box.Y1) { - return IntersectionInside - } - return IntersectionOutSide -} - -func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType { - box := t.BBox() - for _, index := range p.indices { - var intersects bool - index.Search(box, func(item rtree.Item) bool { - ls := item.(lineSegment) - other := make(lineSegment, 4) - for i := range t { - other[0] = t[i].X - other[1] = t[i].Y - other[2] = t[(i+1)%len(t)].X - other[3] = t[(i+1)%len(t)].Y - if ls.intersectsLineSegment(other) { - intersects = true - return false - } - } - return true - }) - if intersects { - return IntersectionOverlaps - } - } - // No intersection -> check inside or outside - // if an abritrary point is inside or not. - pX, pY := t[0].X, t[0].Y - - // Check holes first: inside a hole means outside. - if len(p.rings) > 1 { - for _, hole := range p.rings[1:] { - if contains(hole, pX, pY) { - return IntersectionOutSide - } - } - } - - // Check shell - if contains(p.rings[0], pX, pY) { - return IntersectionInside - } - return IntersectionOutSide -} - -func (rng ring) length() int { - return len(rng) / 2 -} - -func (rng ring) point(i int) (float64, float64) { - i *= 2 - return rng[i], rng[i+1] -} - -type segments interface { - length() int - point(int) (float64, float64) -} - -func contains(s segments, pX, pY float64) bool { - - n := s.length() - if n < 3 { - return false - } - - sX, sY := s.point(0) - eX, eY := s.point(n - 1) - - const eps = 0.0000001 - - if math.Abs(sX-eX) > eps || math.Abs(sY-eY) > eps { - // It's not closed! - return false - } - - var inside bool - - for i := 1; i < n; i++ { - eX, eY := s.point(i) - if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) { - inside = !inside - } - sX, sY = eX, eY - } - - return inside -} - -// 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(pX, pY, sX, sY, eX, eY float64) bool { - - // Always ensure that the the first point - // has a y coordinate that is less than the second point - if sY > eY { - // Switch the points if otherwise. - sX, sY, eX, eY = eX, eY, sX, sY - } - - // Move the point's y coordinate - // outside of the bounds of the testing region - // so we can start drawing a ray - for pY == sY || pY == eY { - pY = math.Nextafter(pY, math.Inf(1)) - } - - // If we are outside of the polygon, indicate so. - if pY < sY || pY > eY { - return false - } - - if sX > eX { - if pX > sX { - return false - } - if pX < eX { - return true - } - } else { - if pX > eX { - return false - } - if pX < sX { - return true - } - } - - raySlope := (pY - sY) / (pX - sX) - diagSlope := (eY - sY) / (eX - sX) - - return raySlope >= diagSlope -} - -func (p *Polygon) NumVertices(ring int) int { - if ring < 0 || ring >= len(p.rings) { - return 0 - } - return len(p.rings[ring]) / 2 -} - -func (p *Polygon) Vertices(ring int, fn func(float64, float64)) { - if ring < 0 || ring >= len(p.rings) { - return - } - rng := p.rings[ring] - for i := 0; i < len(rng); i += 2 { - fn(rng[i+0], rng[i+1]) - } -} - -func (p *Polygon) AsWKB() []byte { - - size := 1 + 4 + 4 - for _, r := range p.rings { - size += 4 + len(r)*8 - } - - buf := bytes.NewBuffer(make([]byte, 0, size)) - - binary.Write(buf, binary.LittleEndian, wkb.NDR) - binary.Write(buf, binary.LittleEndian, wkb.Polygon) - binary.Write(buf, binary.LittleEndian, uint32(len(p.rings))) - - for _, r := range p.rings { - binary.Write(buf, binary.LittleEndian, uint32(len(r)/2)) - for i := 0; i < len(r); i += 2 { - binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0])) - binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1])) - } - } - - return buf.Bytes() -} - -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) - - log.Printf("info: Number of rings: %d\n", len(rngs)) - - for rng := uint32(0); rng < numRings; rng++ { - var numVertices uint32 - if err = binary.Read(r, order, &numVertices); err != nil { - return err - } - - log.Printf("info: Number of vertices in ring %d: %d\n", rng, numVertices) - - 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 -} - -func (p *Polygon) FromLineStringZ(ls LineStringZ) { - r := make([]float64, 2*len(ls)) - var pos int - for i := range ls { - r[pos+0] = ls[i].X - r[pos+1] = ls[i].Y - pos += 2 - } - p.rings = []ring{r} -}
--- a/pkg/octree/raster.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,404 +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) 2019 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 ( - "log" - "math" - "runtime" - "sync" - "time" - - "gemma.intevation.de/gemma/pkg/common" - "gemma.intevation.de/gemma/pkg/wkb" - "github.com/fogleman/contourmap" -) - -type Raster struct { - BBox Box2D - CellSize float64 - XCells int - YCells int - Cells []float64 -} - -const noData = -math.MaxFloat64 - -func NewRaster(bbox Box2D, cellSize float64) *Raster { - - width, height := bbox.Size() - - log.Printf("info raster extent: %.2f / %.2f", width, height) - - xCells := int(math.Ceil(width / cellSize)) - yCells := int(math.Ceil(height / cellSize)) - - log.Printf("info raster size: %d / %d\n", xCells, yCells) - - size := (xCells + 2) * (yCells + 2) - cells := make([]float64, size) - for i := range cells { - cells[i] = noData - } - return &Raster{ - BBox: bbox, - CellSize: cellSize, - XCells: xCells, - YCells: yCells, - Cells: cells, - } -} - -func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) { - var wg sync.WaitGroup - - rows := make(chan int) - - rasterRow := func() { - defer wg.Done() - quat := 0.25 * r.CellSize - for i := range rows { - pos := (i+1)*(r.XCells+2) + 1 - row := r.Cells[pos : pos+r.XCells] - py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 - px := r.BBox.X1 + r.CellSize/2 - y1 := py - quat - y2 := py + quat - for j := range row { - var n int - var sum float64 - - if v, ok := eval(px-quat, y1); ok { - sum = v - n = 1 - } - if v, ok := eval(px-quat, y2); ok { - sum += v - n++ - } - if v, ok := eval(px+quat, y1); ok { - sum += v - n++ - } - if v, ok := eval(px+quat, y2); ok { - sum += v - n++ - } - - if n > 0 { - row[j] = sum / float64(n) - } - px += r.CellSize - } - } - } - - for n := runtime.NumCPU(); n >= 1; n-- { - wg.Add(1) - go rasterRow() - } - - for i := 0; i < r.YCells; i++ { - rows <- i - } - close(rows) -} - -func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) { - var wg sync.WaitGroup - - rows := make(chan int) - - rasterRow := func() { - defer wg.Done() - quat := 0.25 * r.CellSize - for i := range rows { - pos := (i+1)*(r.XCells+2) + 1 - row := r.Cells[pos : pos+r.XCells] - py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 - px := r.BBox.X1 + r.CellSize/2 - y1 := py - quat - y2 := py + quat - for j, old := range row { - // only diff where values - if old == noData { - px += r.CellSize - continue - } - var n int - var sum float64 - - if v, ok := eval(px-quat, y1); ok { - sum = v - n = 1 - } - if v, ok := eval(px-quat, y2); ok { - sum += v - n++ - } - if v, ok := eval(px+quat, y1); ok { - sum += v - n++ - } - if v, ok := eval(px+quat, y2); ok { - sum += v - n++ - } - - if n > 0 { - row[j] -= sum / float64(n) - } else { - row[j] = noData - } - - px += r.CellSize - } - } - } - - for n := runtime.NumCPU(); n >= 1; n-- { - wg.Add(1) - go rasterRow() - } - - for i := 0; i < r.YCells; i++ { - rows <- i - } - close(rows) -} - -func (r *Raster) ZExtent() (float64, float64, bool) { - min, max := math.MaxFloat64, -math.MaxFloat64 - for _, v := range r.Cells { - if v == noData { - continue - } - if v < min { - min = v - } - if v > max { - max = v - } - } - return min, max, min != math.MaxFloat64 -} - -func (r *Raster) Trace(heights []float64) []wkb.MultiPolygonGeom { - start := time.Now() - - tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells) - - areas := make([]wkb.MultiPolygonGeom, len(heights)) - - reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize) - reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize) - - var wg sync.WaitGroup - - cnts := make(chan int) - - doContours := func() { - defer wg.Done() - for hIdx := range cnts { - if c := tracer.Contours(heights[hIdx]); len(c) > 0 { - areas[hIdx] = buildMultipolygon(c, reprojX, reprojY) - } - } - } - - for n := runtime.NumCPU(); n >= 1; n-- { - wg.Add(1) - go doContours() - } - - for i := range heights { - cnts <- i - } - close(cnts) - - wg.Wait() - log.Printf("info: Tracing areas took %v\n", time.Since(start)) - - return areas -} - -type contour []contourmap.Point - -type bboxNode struct { - box Box2D - cnt contour - children []*bboxNode -} - -func (cnt contour) contains(o contour) bool { - return contains(cnt, o[0].X, o[0].Y) || - contains(cnt, o[len(o)/2].X, o[len(o)/2].Y) -} - -func (cnt contour) length() int { - return len(cnt) -} - -func (cnt contour) point(i int) (float64, float64) { - return cnt[i].X, cnt[i].Y -} - -func (cnt contour) bbox() Box2D { - minX, minY := math.MaxFloat64, math.MaxFloat64 - maxX, maxY := -math.MaxFloat64, -math.MaxFloat64 - - for _, p := range cnt { - if p.X < minX { - minX = p.X - } - if p.X > maxX { - maxX = p.X - } - if p.Y < minY { - minY = p.Y - } - if p.Y > maxY { - maxY = p.Y - } - } - return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY} -} - -func (bn *bboxNode) insert(cnt contour, box Box2D) { - // check if children are inside new - var nr *bboxNode - - for i, r := range bn.children { - if r.box.Inside(box) && cnt.contains(r.cnt) { - if nr == nil { - nr = &bboxNode{box: box, cnt: cnt} - } - nr.children = append(nr.children, r) - bn.children[i] = nil - } - } - - // we have a new child - if nr != nil { - // compact the list - for i := len(bn.children) - 1; i >= 0; i-- { - if bn.children[i] == nil { - if i < len(bn.children)-1 { - copy(bn.children[i:], bn.children[i+1:]) - } - bn.children[len(bn.children)-1] = nil - bn.children = bn.children[:len(bn.children)-1] - } - } - bn.children = append(bn.children, nr) - return - } - - // check if new is inside an old - for _, r := range bn.children { - if box.Inside(r.box) && r.cnt.contains(cnt) { - r.insert(cnt, box) - return - } - } - - // its a new child node. - nr = &bboxNode{box: box, cnt: cnt} - bn.children = append(bn.children, nr) -} - -func (bn *bboxNode) insertRoot(cnt contour) { - bn.insert(cnt, cnt.bbox()) -} - -type bboxOutFunc func(contour, []contour) - -func (bn *bboxNode) generate(out bboxOutFunc) { - - var grands []*bboxNode - - holes := make([]contour, len(bn.children)) - - for i, ch := range bn.children { - holes[i] = ch.cnt - grands = append(grands, ch.children...) - } - out(bn.cnt, holes) - - // the grand children are new polygons. - for _, grand := range grands { - grand.generate(out) - } -} - -func (bn *bboxNode) generateRoot(out bboxOutFunc) { - for _, r := range bn.children { - r.generate(out) - } -} - -func buildMultipolygon( - cnts []contourmap.Contour, - reprojX, reprojY func(float64) float64, -) wkb.MultiPolygonGeom { - - var forest bboxNode - - for _, cnt := range cnts { - forest.insertRoot(contour(cnt)) - } - - //log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots)) - - var mp wkb.MultiPolygonGeom - - out := func(sh contour, hls []contour) { - - polygon := make(wkb.PolygonGeom, 1+len(hls)) - - // Handle shell - shell := make(wkb.LinearRingGeom, len(sh)) - for j, pt := range sh { - shell[j] = wkb.PointGeom{ - X: reprojX(pt.X), - Y: reprojY(pt.Y), - } - } - if shell.CCW() { - shell.Reverse() - } - polygon[0] = shell - - // handle holes - for i, hl := range hls { - hole := make(wkb.LinearRingGeom, len(hl)) - for j, pt := range hl { - hole[j] = wkb.PointGeom{ - X: reprojX(pt.X), - Y: reprojY(pt.Y), - } - } - if !hole.CCW() { - hole.Reverse() - } - polygon[1+i] = hole - } - - mp = append(mp, polygon) - } - - forest.generateRoot(out) - - return mp -}
--- a/pkg/octree/strtree.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,486 +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 - -import ( - "bytes" - "compress/gzip" - "encoding/binary" - "io" - "log" - "math" - "sort" -) - -const STRTreeDefaultEntries = 8 - -type STRTree struct { - Entries int - tin *Tin - index []int32 - 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 - - for len(stack) > 0 { - top := stack[len(stack)-1] - stack = stack[:len(stack)-1] - - if top > 0 { // node - if s.bbox(top).Contains(x, y) { - n := s.index[top+1] - stack = append(stack, s.index[top+2:top+2+n]...) - } - } else { // leaf - top = -top - 1 - if !s.bbox(top).Contains(x, y) { - 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]], - } - if t.Contains(x, y) { - return t.Plane3D().Z(x, y), true - } - } - } - } - return 0, false -} - -func (s *STRTree) Build(t *Tin) { - - if s.Entries == 0 { - s.Entries = STRTreeDefaultEntries - } - - s.tin = t - - all := make([]int32, len(t.Triangles)) - - for i := range all { - all[i] = int32(i) - } - - s.index = append(s.index, 0) - - root := s.build(all) - - s.index[0] = root -} - -func (s *STRTree) BuildWithout(t *Tin, remove map[int32]struct{}) { - - if s.Entries == 0 { - s.Entries = STRTreeDefaultEntries - } - - s.tin = t - - all := make([]int32, 0, len(t.Triangles)-len(remove)) - - for i := 0; i < len(t.Triangles); i++ { - idx := int32(i) - if _, found := remove[idx]; !found { - all = append(all, idx) - } - } - - s.index = append(s.index, 0) - - root := s.build(all) - - s.index[0] = root -} - -func (s *STRTree) Clip(p *Polygon) map[int32]struct{} { - - removed := make(map[int32]struct{}) - - stack := []int32{s.index[0]} - - vertices := s.tin.Vertices - - for len(stack) > 0 { - top := stack[len(stack)-1] - stack = stack[:len(stack)-1] - - if top > 0 { // node - switch p.IntersectionBox2D(s.bbox(top)) { - case IntersectionInside: - // all triangles are inside polygon - case IntersectionOutSide: - // all triangles are outside polygon - s.allTriangles(top, removed) - default: - // mixed bag - n := s.index[top+1] - stack = append(stack, s.index[top+2:top+2+n]...) - } - } else { // leaf - top = -top - 1 - switch p.IntersectionBox2D(s.bbox(top)) { - case IntersectionInside: - // all triangles are inside polygon - case IntersectionOutSide: - // all triangles are outside polygon - n := s.index[top+1] - for _, idx := range s.index[top+2 : top+2+n] { - removed[idx] = struct{}{} - } - default: - 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]], - } - if p.IntersectionWithTriangle(&t) != IntersectionInside { - removed[idx] = struct{}{} - } - } - } - } - } - - 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 { - top := stack[len(stack)-1] - stack = stack[:len(stack)-1] - if top > 0 { // node - n := s.index[top+1] - stack = append(stack, s.index[top+2:top+2+n]...) - } else { // leaf - top = -top - 1 - n := s.index[top+1] - for _, idx := range s.index[top+2 : top+2+n] { - tris[idx] = struct{}{} - } - } - } -} - -func (s *STRTree) build(items []int32) int32 { - sort.Slice(items, func(i, j int) bool { - ti := s.tin.Triangles[items[i]] - tj := s.tin.Triangles[items[j]] - return s.tin.Vertices[ti[0]].X < s.tin.Vertices[tj[0]].X - }) - - P := int(math.Ceil(float64(len(items)) / float64(s.Entries))) - S := int(math.Ceil(math.Sqrt(float64(P)))) - - slices := s.strSplit(items, S) - - leaves := s.strJoin( - slices, S, - func(i, j int32) bool { - ti := s.tin.Triangles[i] - tj := s.tin.Triangles[j] - return s.tin.Vertices[ti[0]].Y < s.tin.Vertices[tj[0]].Y - }, - s.allocLeaf, - ) - - return s.buildNodes(leaves) -} - -func (s *STRTree) buildNodes(items []int32) int32 { - - if len(items) <= s.Entries { - return s.allocNode(items) - } - - sort.Slice(items, func(i, j int) bool { - return s.bbox(items[i]).X1 < s.bbox(items[j]).X1 - }) - - P := int(math.Ceil(float64(len(items)) / float64(s.Entries))) - S := int(math.Ceil(math.Sqrt(float64(P)))) - - slices := s.strSplit(items, S) - - nodes := s.strJoin( - slices, S, - func(i, j int32) bool { return s.bbox(i).Y1 < s.bbox(j).Y1 }, - s.allocNode, - ) - - return s.buildNodes(nodes) -} - -func (s *STRTree) bbox(idx int32) Box2D { - if idx < 0 { // Don't care if leaf or node. - idx = -idx - 1 - } - return s.bboxes[s.index[idx]] -} - -func (s *STRTree) strSplit(items []int32, S int) [][]int32 { - sm := S * s.Entries - slices := make([][]int32, S) - for i := range slices { - var n int - if l := len(items); l < sm { - n = l - } else { - n = sm - } - slices[i] = items[:n] - items = items[n:] - } - return slices -} - -func (s *STRTree) strJoin( - slices [][]int32, S int, - less func(int32, int32) bool, - alloc func([]int32) int32, -) []int32 { - nodes := make([]int32, 0, S*S) - - for _, slice := range slices { - sort.Slice(slice, func(i, j int) bool { - return less(slice[i], slice[j]) - }) - - for len(slice) > 0 { - var n int - if l := len(slice); l >= s.Entries { - n = s.Entries - } else { - n = l - } - nodes = append(nodes, alloc(slice[:n])) - slice = slice[n:] - } - } - return nodes -} - -func (s *STRTree) allocNode(items []int32) int32 { - pos := len(s.index) - s.index = append(s.index, 0, int32(len(items))) - s.index = append(s.index, items...) - if len(items) > 0 { - box := s.bbox(items[0]) - for i := 1; i < len(items); i++ { - box = box.Union(s.bbox(items[i])) - } - s.index[pos] = int32(s.allocBBox(box)) - } - return int32(pos) -} - -func (s *STRTree) allocBBox(box Box2D) int { - pos := len(s.bboxes) - s.bboxes = append(s.bboxes, box) - return pos -} - -func (s *STRTree) allocLeaf(items []int32) int32 { - pos := len(s.index) - s.index = append(s.index, 0, int32(len(items))) - s.index = append(s.index, items...) - if len(items) > 0 { - vertices := s.tin.Vertices - ti := s.tin.Triangles[items[0]] - t := Triangle{ - vertices[ti[0]], - vertices[ti[1]], - vertices[ti[2]], - } - box := t.BBox() - for i := 1; i < len(items); i++ { - it := items[i] - ti := s.tin.Triangles[it] - t := Triangle{ - vertices[ti[0]], - vertices[ti[1]], - vertices[ti[2]], - } - box = box.Union(t.BBox()) - } - s.index[pos] = int32(s.allocBBox(box)) - } - return int32(-(pos + 1)) -}
--- a/pkg/octree/tin.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,271 +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 - -import ( - "bytes" - "encoding/binary" - "errors" - "fmt" - "io" - "log" - "math" - - "gemma.intevation.de/gemma/pkg/models" - "gemma.intevation.de/gemma/pkg/wkb" -) - -var ( - errNoByteSlice = errors.New("Not a byte slice") - errTooLessPoints = errors.New("Too less points") -) - -// Tin stores a mesh of triangles with common vertices. -type Tin struct { - // EPSG holds the projection. - EPSG uint32 - // Vertices are the shared vertices. - Vertices []Vertex - // Triangles are the triangles. - Triangles [][]int32 - - // Min is the lower left corner of the bbox. - Min Vertex - // Max is the upper right corner of the bbox. - Max Vertex -} - -func (t *Tin) Clip(polygon *Polygon) map[int32]struct{} { - var tree STRTree - tree.Build(t) - return tree.Clip(polygon) -} - -// FromWKB constructs the TIN from a WKB representation. -// Shared vertices are identified and referenced by the -// same index. -func (t *Tin) FromWKB(data []byte) error { - log.Printf("info: data length %d\n", len(data)) - - 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.TinZ: - return fmt.Errorf("unknown geometry type %x", geomType) - } - - var num uint32 - if err = binary.Read(r, order, &num); err != nil { - return err - } - - vertices := make([]Vertex, 0, 100000) - - var v Vertex - - v2i := make(map[Vertex]int32, 100000) - - var indexPool []int32 - - allocIndices := func() []int32 { - if len(indexPool) == 0 { - indexPool = make([]int32, 3*8*1024) - } - ids := indexPool[:3] - indexPool = indexPool[3:] - return ids - } - - var triangles [][]int32 - - min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} - max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} - - for i := uint32(0); i < num; i++ { - - endian, err = r.ReadByte() - 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) - } - - err = binary.Read(r, order, &geomType) - switch { - case err != nil: - return err - case geomType != wkb.TriangleZ: - return fmt.Errorf("unknown geometry type %d", geomType) - } - - var rings uint32 - if err = binary.Read(r, order, &rings); err != nil { - return err - } - triangle := allocIndices() - - for ring := uint32(0); ring < rings; ring++ { - var npoints uint32 - if err = binary.Read(r, order, &npoints); err != nil { - return err - } - - if npoints < 3 { - return errTooLessPoints - } - - for p := uint32(0); p < npoints; p++ { - var x, y, z uint64 - for _, addr := range []*uint64{&x, &y, &z} { - if err = binary.Read(r, order, addr); err != nil { - return err - } - } - if p >= 3 || ring >= 1 { - // Don't store the forth point. - continue - } - // Do this conversion later to spare reflect calls - // and allocs in binary.Read. - v.X = math.Float64frombits(x) - v.Y = math.Float64frombits(y) - v.Z = math.Float64frombits(z) - idx, found := v2i[v] - if !found { - idx = int32(len(vertices)) - v2i[v] = idx - vertices = append(vertices, v) - min.Minimize(v) - max.Maximize(v) - } - triangle[p] = idx - } - } - triangles = append(triangles, triangle) - } - - log.Printf("info: bbox: [[%f, %f], [%f, %f]]\n", - min.X, min.Y, max.X, max.Y) - - *t = Tin{ - EPSG: models.WGS84, - Vertices: vertices, - Triangles: triangles, - Min: min, - Max: max, - } - - return nil -} - -// Scan implements the sql.Scanner interface. -func (t *Tin) Scan(raw interface{}) error { - if raw == nil { - return nil - } - data, ok := raw.([]byte) - if !ok { - return errNoByteSlice - } - return t.FromWKB(data) -} - -func (t *Tin) serialize(w io.Writer) error { - - if err := binary.Write(w, binary.LittleEndian, t.EPSG); err != nil { - return err - } - - if err := t.Min.Write(w); err != nil { - return err - } - if err := t.Max.Write(w); err != nil { - return err - } - - if err := binary.Write( - w, binary.LittleEndian, uint32(len(t.Vertices))); err != nil { - 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 { - 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) - - if err := binary.Write( - w, binary.LittleEndian, uint32(len(t.Triangles))); err != nil { - return err - } - - var buf [binary.MaxVarintLen32]byte - var written int - var last int32 - for _, triangle := range t.Triangles { - for _, idx := range triangle { - value := idx - last - n := binary.PutVarint(buf[:], int64(value)) - 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 = idx - } - } - log.Printf("info: compressed tin indices in bytes: %d (%d)\n", - written, 3*4*len(t.Triangles)) - - return nil -}
--- a/pkg/octree/triangulation.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,318 +0,0 @@ -// 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" - "log" - "math" - - "gonum.org/v1/gonum/stat" -) - -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) EstimateTooLong() float64 { - - num := len(t.Triangles) / 3 - - lengths := make([]float64, 0, num) - - points := t.Points - -tris: - for i := 0; i < num; i++ { - idx := i * 3 - var max float64 - vs := t.Triangles[idx : idx+3] - for j, vj := range vs { - if t.Halfedges[idx+j] < 0 { - continue tris - } - k := (j + 1) % 3 - if l := points[vj].Distance2D(points[vs[k]]); l > max { - max = l - } - } - lengths = append(lengths, max) - } - - std := stat.StdDev(lengths, nil) - return 3.5 * std -} - -func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) { - - if tooLong <= 0 { - tooLong = t.EstimateTooLong() - } - - tooLong *= tooLong - - var candidates []int32 - - points := t.Points - - for i, num := 0, len(t.Triangles)/3; i < num; i++ { - idx := i * 3 - var max float64 - vs := t.Triangles[idx : idx+3] - for j, vj := range vs { - k := (j + 1) % 3 - if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max { - max = l - } - } - if max > tooLong { - candidates = append(candidates, int32(i)) - } - } - - removed := map[int32]struct{}{} - - isBorder := func(n int32) bool { - n *= 3 - for i := int32(0); i < 3; i++ { - e := n + i - o := t.Halfedges[e] - if o < 0 { - return true - } - if _, found := removed[o/3]; found { - return true - } - } - return false - } - - var newCandidates []int32 - - log.Printf("info: candidates: %d\n", len(candidates)) - for len(candidates) > 0 { - - oldRemoved := len(removed) - - for _, i := range candidates { - - if isBorder(i) { - removed[i] = struct{}{} - } else { - newCandidates = append(newCandidates, i) - } - } - - if oldRemoved == len(removed) { - break - } - - candidates = newCandidates - newCandidates = newCandidates[:0] - } - - log.Printf("info: candidates left: %d\n", len(candidates)) - log.Printf("info: triangles: %d\n", len(t.Triangles)/3) - log.Printf("info: triangles to remove: %d\n", len(removed)) - - type edge struct { - a, b int32 - prev, next *edge - } - - isClosed := func(e *edge) bool { - for curr := e.next; curr != nil; curr = curr.next { - if curr == e { - return true - } - } - return false - } - - open := map[int32]*edge{} - var rings []*edge - - for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ { - if _, found := removed[i]; found { - continue - } - n := i * 3 - for j := int32(0); j < 3; j++ { - e := n + j - f := t.Halfedges[e] - if f >= 0 { - if _, found := removed[f/3]; !found { - continue - } - } - a := t.Triangles[e] - b := t.Triangles[n+(j+1)%3] - - curr := &edge{a: a, b: b} - - if old := open[a]; old != nil { - delete(open, a) - if old.a == a { - old.prev = curr - curr.next = old - } else { - old.next = curr - curr.prev = old - } - - if isClosed(old) { - rings = append(rings, old) - } - } else { - open[a] = curr - } - - if old := open[b]; old != nil { - delete(open, b) - if old.b == b { - old.next = curr - curr.prev = old - } else { - old.prev = curr - curr.next = old - } - - if isClosed(old) { - rings = append(rings, old) - } - } else { - open[b] = curr - } - } - } - - if len(open) > 0 { - log.Printf("warn: open vertices left: %d\n", len(open)) - } - - if len(rings) == 0 { - log.Println("warn: no ring found") - return nil, removed - } - - curr := rings[0] - - polygon := LineStringZ{ - points[curr.a], - points[curr.b], - } - - for curr = curr.next; curr != rings[0]; curr = curr.next { - polygon = append(polygon, points[curr.b]) - } - - polygon = append(polygon, t.Points[rings[0].a]) - - log.Printf("length of boundary: %d\n", len(polygon)) - - return polygon, removed -} - -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 -}
--- a/pkg/octree/triangulator.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,375 +0,0 @@ -// 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 ( - "errors" - "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 (tri *triangulator) Len() int { - return len(tri.points) -} - -func (tri *triangulator) Swap(i, j int) { - tri.ids[i], tri.ids[j] = tri.ids[j], tri.ids[i] -} - -func (tri *triangulator) Less(i, j int) bool { - d1 := tri.squaredDistances[tri.ids[i]] - d2 := tri.squaredDistances[tri.ids[j]] - if d1 != d2 { - return d1 < d2 - } - p1 := tri.points[tri.ids[i]] - p2 := tri.points[tri.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 errors.New("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 (tri *triangulator) hashKey(point Vertex) int { - d := point.Sub2D(tri.center) - return int(pseudoAngle(d.X, d.Y) * float64(len(tri.hash))) -} - -func (tri *triangulator) hashEdge(e *node) { - tri.hash[tri.hashKey(tri.points[e.i])] = e -} - -func (tri *triangulator) addTriangle(i0, i1, i2, a, b, c int32) int32 { - i := int32(tri.trianglesLen) - tri.triangles[i] = i0 - tri.triangles[i+1] = i1 - tri.triangles[i+2] = i2 - tri.link(i, a) - tri.link(i+1, b) - tri.link(i+2, c) - tri.trianglesLen += 3 - return i -} - -func (tri *triangulator) link(a, b int32) { - tri.halfedges[a] = b - if b >= 0 { - tri.halfedges[b] = a - } -} - -func (tri *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 := tri.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 := tri.triangles[ar] - pr := tri.triangles[a] - pl := tri.triangles[al] - p1 := tri.triangles[bl] - - illegal := inCircle(tri.points[p0], tri.points[pr], tri.points[pl], tri.points[p1]) - - if illegal { - tri.triangles[a] = p1 - tri.triangles[b] = p0 - - // edge swapped on the other side of the hull (rare) - // fix the halfedge reference - if tri.halfedges[bl] == -1 { - e := tri.hull - for { - if e.t == bl { - e.t = a - break - } - e = e.next - if e == tri.hull { - break - } - } - } - - tri.link(a, tri.halfedges[bl]) - tri.link(b, tri.halfedges[ar]) - tri.link(ar, bl) - - br := b0 + (b+1)%3 - - tri.legalize(a) - return tri.legalize(br) - } - - return ar -} - -func (tri *triangulator) convexHull() []Vertex { - var result []Vertex - e := tri.hull - for e != nil { - result = append(result, tri.points[e.i]) - e = e.prev - if e == tri.hull { - break - } - } - return result -}
--- a/pkg/octree/util.go Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -// 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 Tue Nov 05 13:01:37 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1229 +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 - -import ( - "bytes" - "encoding/binary" - "fmt" - "io" - "log" - "math" - "sort" - - "gemma.intevation.de/gemma/pkg/wkb" -) - -type ( - Point struct { - X float64 - Y float64 - } - - // Vertex is a 3D vertex. - Vertex struct { - X float64 - Y float64 - Z float64 - } - - // Triangle is a triangle consisting of three vertices. - Triangle [3]Vertex - - // Line is a line defined by first vertex on that line - // and the second being the direction. - Line [2]Vertex - - // Box is a 3D box. - Box [2]Vertex - - // MultiPointZ is a set of vertices. - MultiPointZ []Vertex - - // LineStringZ is a line string formed of vertices. - LineStringZ []Vertex - - // MultiLineStringZ is a set of line strings. - MultiLineStringZ []LineStringZ - - // Box2D is 2D area from (X1, Y1) to (X2, Y2). - Box2D struct { - X1 float64 - Y1 float64 - X2 float64 - Y2 float64 - } - - // Plane2D is a 2D plane (a line in 2D space). - Plane2D struct { - A float64 - 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 (t *Triangle) BBox() Box2D { - minX := math.Min(math.Min(t[0].X, t[1].X), t[2].X) - maxX := math.Max(math.Max(t[0].X, t[1].X), t[2].X) - minY := math.Min(math.Min(t[0].Y, t[1].Y), t[2].Y) - maxY := math.Max(math.Max(t[0].Y, t[1].Y), t[2].Y) - return Box2D{ - X1: minX, Y1: minY, - X2: maxX, Y2: maxY, - } -} - -func (a Box2D) Inside(b Box2D) bool { - return a.X1 >= b.X1 && a.X2 <= b.X2 && - a.Y1 >= b.Y1 && a.Y2 <= b.Y2 -} - -func (a Box2D) Size() (float64, float64) { - return a.X2 - a.X1, a.Y2 - a.Y1 -} - -func (a Box2D) Empty() bool { - const eps = 0.0000001 - return math.Abs(a.X2-a.X1) < eps && - math.Abs(a.Y2-a.Y1) < eps -} - -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 (p Plane3D) Eval(v Vertex) float64 { - return p.A*v.X + p.B*v.Y + p.C*v.Z + p.D -} - -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) -} - -func (v Vertex) Distance(w Vertex) float64 { - v = v.Sub(w) - 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) { - if w.X < v.X { - v.X = w.X - } - if w.Y < v.Y { - v.Y = w.Y - } - if w.Z < v.Z { - v.Z = w.Z - } -} - -// Maximize adjust this vertex v to hold the maximum -// values component-wise of v and w. -func (v *Vertex) Maximize(w Vertex) { - if w.X > v.X { - v.X = w.X - } - if w.Y > v.Y { - v.Y = w.Y - } - if w.Z > v.Z { - v.Z = w.Z - } -} - -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{ - v.X - w.X, - v.Y - w.Y, - v.Z - w.Z, - } -} - -// Add returns (v + w) component-wise. -func (v Vertex) Add(w Vertex) Vertex { - return Vertex{ - v.X + w.X, - v.Y + w.Y, - v.Z + w.Z, - } -} - -// Scale returns s*v component-wise. -func (v Vertex) Scale(s float64) Vertex { - return Vertex{ - s * v.X, - s * v.Y, - s * v.Z, - } -} - -// Interpolate returns a function that return s*b[1] + b[0] -// component-wise. -func (b Box) Interpolate() func(Vertex) Vertex { - v1, v2 := b[0], b[1] - v2 = v2.Sub(v1) - return func(s Vertex) Vertex { - return Vertex{ - v2.X*s.X + v1.X, - v2.Y*s.Y + v1.Y, - v2.Z*s.Z + v1.Z, - } - } -} - -func (b Box) HasX() bool { return math.Abs(b[0].X-b[1].X) > epsPlane } -func (b Box) HasY() bool { return math.Abs(b[0].Y-b[1].Y) > epsPlane } -func (b Box) HasZ() bool { return math.Abs(b[0].Z-b[1].Z) > epsPlane } - -// Less returns if one of v component is less than the -// corresponing component in w. -func (v Vertex) Less(w Vertex) bool { - return v.X < w.X || v.Y < w.Y || v.Z < w.Z -} - -// NewLine return a line of point/direction. -func NewLine(p1, p2 Vertex) Line { - return Line{ - p2.Sub(p1), - p1, - } -} - -// Eval returns the vertex for t*l[0] + l[1]. -func (l Line) Eval(t float64) Vertex { - return l[0].Scale(t).Add(l[1]) -} - -// IntersectHorizontal returns the intersection point -// for a given z value. -func (l Line) IntersectHorizontal(h float64) Vertex { - t := (h - l[1].Z) / l[0].Z - return l.Eval(t) -} - -func side(z, h float64) int { - switch { - case z < h: - return -1 - case z > h: - return +1 - } - 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].Sub2D(t[0]) - v1 := t[1].Sub2D(t[0]) - v2 := Vertex{X: x, Y: y}.Sub2D(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), -// one vertex (only touching an vertex) or -// two vertices (real intersection). -func (t *Triangle) IntersectHorizontal(h float64) LineStringZ { - sides := [3]int{ - side(t[0].Z, h), - side(t[1].Z, h), - side(t[2].Z, h), - } - - var points LineStringZ - - for i := 0; i < 3; i++ { - j := (i + 1) % 3 - si := sides[i] - sj := sides[j] - - switch { - case si == 0 && sj == 0: - // both on plane - points = append(points, t[i], t[j]) - case si == 0 && sj != 0: - // first on plane - points = append(points, t[i]) - case si != 0 && sj == 0: - // second on plane - points = append(points, t[j]) - case si == sj: - // both on same side - default: - // real intersection - v := NewLine(t[i], t[j]).IntersectHorizontal(h) - points = append(points, v) - } - } - - return points -} - -func linearScale(x1, y1, x2, y2 float64) func(Vertex) float64 { - dx := x2 - x1 - dy := y2 - y1 - - switch { - case dx != 0: - return func(v Vertex) float64 { - return (v.X - x1) / dx - } - case dy != 0: - return func(v Vertex) float64 { - return (v.Y - y1) / dy - } - default: - return func(Vertex) float64 { - return 0 - } - } -} - -func (ls LineStringZ) BBox() Box2D { - - min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} - max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} - - for _, v := range ls { - min.Minimize(v) - max.Maximize(v) - } - - return Box2D{ - X1: min.X, - Y1: min.Y, - X2: max.X, - Y2: max.Y, - } -} - -func (ls LineStringZ) Area() float64 { - return polygonArea(ls) -} - -func (ls LineStringZ) Reverse() { - for i, j := 0, len(ls)-1; i < j; i, j = i+1, j-1 { - ls[i], ls[j] = ls[j], ls[i] - } -} - -func (ls LineStringZ) order(position func(Vertex) float64) { - type posVertex struct { - pos float64 - v Vertex - } - positions := make([]posVertex, len(ls)) - for i, v := range ls { - positions[i] = posVertex{position(v), v} - } - sort.Slice(positions, func(i, j int) bool { - return positions[i].pos < positions[j].pos - }) - for i := range positions { - ls[i] = positions[i].v - } -} - -// EpsEquals returns true if v and w are equal component-wise -// with the values within a epsilon range. -func (v Vertex) EpsEquals(w Vertex) bool { - const eps = 1e-5 - return math.Abs(v.X-w.X) < eps && - math.Abs(v.Y-w.Y) < eps && math.Abs(v.Z-w.Z) < eps -} - -// EpsEquals returns true if v and w are equal component-wise -// in the X/Y plane with the values within a epsilon range. -func (v Vertex) EpsEquals2D(w Vertex) bool { - const eps = 1e-5 - return math.Abs(v.X-w.X) < eps && - math.Abs(v.Y-w.Y) < eps -} - -// JoinOnLine joins the the elements of a given multi line string -// under the assumption that the segments are all on the line segment -// from (x1, y1) to (x2, y2). -func (mls MultiLineStringZ) JoinOnLine(x1, y1, x2, y2 float64) MultiLineStringZ { - - position := linearScale(x1, y1, x2, y2) - - type posLineString struct { - pos float64 - line LineStringZ - } - - positions := make([]posLineString, 0, len(mls)) - - for _, ls := range mls { - if len(ls) == 0 { - continue - } - // order the atoms first - ls.order(position) - positions = append(positions, posLineString{position(ls[0]), ls}) - } - - sort.Slice(positions, func(i, j int) bool { - return positions[i].pos < positions[j].pos - }) - - out := make(MultiLineStringZ, 0, len(positions)) - - var ignored int - - for i := range positions { - curr := positions[i].line - if l := len(out); l > 0 { - last := out[l-1] - - if last[len(last)-1].EpsEquals(curr[0]) { - out[l-1] = append(last[:len(last)-1], curr...) - continue - } - if position(last[len(last)-1]) > position(curr[0]) { - ignored++ - continue - } - } - out = append(out, curr) - } - - log.Printf("info: ignored parts: %d\n", ignored) - - return out -} - -// Write writes a Vertex as three 64 bit values in little endian order -// to the given writer. -func (v *Vertex) Write(w io.Writer) error { - if err := binary.Write( - w, binary.LittleEndian, math.Float64bits(v.X)); err != nil { - return err - } - if err := binary.Write( - w, binary.LittleEndian, math.Float64bits(v.Y)); err != nil { - return err - } - return binary.Write( - w, binary.LittleEndian, math.Float64bits(v.Z)) -} - -// Read fills this vertex with three 64 bit values stored as -// little endian from the given reader. -func (v *Vertex) Read(r io.Reader) error { - var buf [8]byte - b := buf[:] - if _, err := io.ReadFull(r, b); err != nil { - return nil - } - v.X = math.Float64frombits(binary.LittleEndian.Uint64(b)) - if _, err := io.ReadFull(r, b); err != nil { - return nil - } - v.Y = math.Float64frombits(binary.LittleEndian.Uint64(b)) - if _, err := io.ReadFull(r, b); err != nil { - return nil - } - v.Z = math.Float64frombits(binary.LittleEndian.Uint64(b)) - return nil -} - -// AsWKB returns the WKB representation of the given multi line string. -func (mls MultiLineStringZ) AsWKB() []byte { - - // pre-calculate size to avoid reallocations. - size := 1 + 4 + 4 - for _, ml := range mls { - size += 1 + 4 + 4 + len(ml)*3*8 - } - - buf := bytes.NewBuffer(make([]byte, 0, size)) - - 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, 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)) - binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) - binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z)) - } - } - - return buf.Bytes() -} - -// AsWKB2D returns the WKB representation of the given multi line string -// leaving the z component out. -func (mls MultiLineStringZ) AsWKB2D() []byte { - - // pre-calculate size to avoid reallocations. - size := 1 + 4 + 4 - for _, ml := range mls { - size += 1 + 4 + 4 + len(ml)*2*8 - } - - buf := bytes.NewBuffer(make([]byte, 0, size)) - - 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, 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)) - binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) - } - } - - return buf.Bytes() -} - -func (ls LineStringZ) CCW() bool { - var sum float64 - for i, v1 := range ls { - v2 := ls[(i+1)%len(ls)] - sum += (v2.X - v1.X) * (v2.Y + v1.Y) - } - return sum > 0 -} - -// Join joins two lines leaving the first of the second out. -func (ls LineStringZ) Join(other LineStringZ) LineStringZ { - nline := make(LineStringZ, len(ls)+len(other)-1) - copy(nline, ls) - copy(nline[len(ls):], other[1:]) - return nline -} - -// Merge merges line segments of a given multi line string -// by finding common start and end vertices. -func (mls MultiLineStringZ) Merge() MultiLineStringZ { - - var out MultiLineStringZ - - min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} - - for _, line := range mls { - for _, v := range line { - min.Minimize(v) - } - } - - type point struct { - x int64 - y int64 - } - - const precision = 1e7 - - quant := func(v Vertex) point { - return point{ - x: int64(math.Round((v.X - min.X) * precision)), - y: int64(math.Round((v.Y - min.Y) * precision)), - } - } - - heads := make(map[point]*[]LineStringZ) - - for _, line := range mls { - if len(line) < 2 { - out = append(out, line) - continue - } - head := quant(line[0]) - tail := quant(line[len(line)-1]) - if head == tail { // its already a ring - out = append(out, line) - continue - } - - if hs := heads[tail]; hs != nil { - l := len(*hs) - last := (*hs)[l-1] - if l == 1 { - delete(heads, tail) - } else { - (*hs)[l-1] = nil - *hs = (*hs)[:l-1] - } - line = line.Join(last) - - if head == quant(line[len(line)-1]) { // its a ring - out = append(out, line) - continue - } - } - - if hs := heads[head]; hs != nil { - *hs = append(*hs, line) - } else { - heads[head] = &[]LineStringZ{line} - } - } - -again: - for head, lines := range heads { - for i, line := range *lines { - tail := quant(line[len(line)-1]) - for hs := heads[tail]; hs != nil && len(*hs) > 0; hs = heads[tail] { - l := len(*hs) - last := (*hs)[l-1] - (*hs)[l-1] = nil - *hs = (*hs)[:l-1] - line = line.Join(last) - - if tail = quant(line[len(line)-1]); head == tail { // its a ring - out = append(out, line) - // remove from current lines - copy((*lines)[i:], (*lines)[i+1:]) - (*lines)[len(*lines)-1] = nil - *lines = (*lines)[:len(*lines)-1] - goto again - } - // overwrite in current lines - (*lines)[i] = line - } - } - } - - // rings := len(out) - - // The rest are open line strings. - for _, lines := range heads { - for _, line := range *lines { - out = append(out, line) - } - } - - // log.Printf("segments before/after merge: %d/%d (%d rings)\n", - // len(mls), len(out), rings) - - 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 { - return a.X1 - } - return a.X2 -} - -// Yi returns the i-th y component. -func (a Box2D) Yi(i int) float64 { - if i == 0 { - return a.Y1 - } - return a.Y2 -} - -func (a Box2D) Union(b Box2D) Box2D { - return Box2D{ - X1: math.Min(a.X1, b.X1), - Y1: math.Min(a.Y1, b.Y1), - X2: math.Max(a.X2, b.X2), - Y2: math.Max(a.Y2, b.Y2), - } -} - -func (a Box2D) Area() float64 { - return (a.X2 - a.X1) * (a.Y2 - a.Y1) -} - -// NewPlane2D creates a new Plane2D from two given points. -func NewPlane2D(x1, y1, x2, y2 float64) Plane2D { - b := x2 - x1 - a := -(y2 - y1) - - l := math.Sqrt(a*a + b*b) - a /= l - b /= l - - // a*x1 + b*y1 + c = 0 - c := -(a*x1 + b*y1) - return Plane2D{a, b, c} -} - -// Eval determines the distance of a given point -// from the plane. The sign of the result indicates -// the sideness. -func (p Plane2D) Eval(x, y float64) float64 { - return p.A*x + p.B*y + p.C -} - -const epsPlane = 1e-5 - -func sides(s int, x float64) int { - if math.Signbit(x) { - return s | 2 - } - return s | 1 -} - -// IntersectsPlane checks if a Box2D intersects with -// a given Plane2D. -func (a Box2D) IntersectsPlane(p Plane2D) bool { - var s int - for i := 0; i < 2; i++ { - x := a.Xi(i) - for j := 0; j < 2; j++ { - y := a.Yi(j) - v := p.Eval(x, y) - if math.Abs(v) < epsPlane { - //log.Printf("on line") - return true - } - if s = sides(s, v); s == 3 { - //log.Printf("... on both sides (djt)") - return true - } - } - } - //log.Printf("side: %d\n", s) - return false -} - -// Cross calculates the cross product of two vertices. -func (v Vertex) Cross(w Vertex) Vertex { - return Vertex{ - v.Y*w.Z - v.Z*w.Y, - v.Z*w.X - v.X*w.Z, - v.X*w.Y - v.Y*w.X, - } -} - -// Intersection calcultes the 2D intersection point of -// two Plane2Ds. If they do not intersect the returned -// bool flags is set to false. -func (p Plane2D) Intersection(o Plane2D) (float64, float64, bool) { - - u1 := Vertex{p.A, p.B, p.C} - u2 := Vertex{o.A, o.B, o.C} - - plane := u1.Cross(u2) - - if plane.Z == 0 { - return 0, 0, false - } - - return plane.X / plane.Z, plane.Y / plane.Z, true -} - -// VerticalLine is a 2D line segment. -type VerticalLine struct { - x1 float64 - y1 float64 - x2 float64 - y2 float64 - - line Plane2D -} - -// NewVerticalLine creates a new 2D line segment. -func NewVerticalLine(x1, y1, x2, y2 float64) *VerticalLine { - return &VerticalLine{ - x1: x1, - y1: y1, - x2: x2, - y2: y2, - line: NewPlane2D(x1, y1, x2, y2), - } -} - -func onPlane(x float64) bool { return math.Abs(x) < epsPlane } - -func relative(v1, v2 Vertex) func(x, y float64) float64 { - ls := linearScale(v1.X, v1.Y, v2.X, v2.Y) - return func(x, y float64) float64 { - return ls(Vertex{x, y, 0}) - } -} - -func interpolate(a, b float64) func(float64) float64 { - return func(x float64) float64 { - return (b-a)*x + a - } -} - -func linear(v1, v2 Vertex) func(t float64) Vertex { - return func(t float64) Vertex { - return Vertex{ - (v2.X-v1.X)*t + v1.X, - (v2.Y-v1.Y)*t + v1.Y, - (v2.Z-v1.Z)*t + v1.Z, - } - } -} - -func inRange(a float64) bool { return 0 <= a && a <= 1 } - -// Intersection intersects a line segment with a triangle. -func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ { - - var out LineStringZ - - //defer func() { log.Printf("length out: %d\n", len(out)) }() - -edges: - for i := 0; i < 3 && len(out) < 2; i++ { - j := (i + 1) % 3 - edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y) - - s1 := edge.Eval(vl.x1, vl.y1) - s2 := edge.Eval(vl.x2, vl.y2) - - o1 := onPlane(s1) - o2 := onPlane(s2) - - // log.Printf("s1, s2: %t %t (%f %f)\n", o1, o2, s1, s2) - - switch { - case o1 && o2: - pos := relative(t[i], t[j]) - t1 := pos(vl.x1, vl.y1) - t2 := pos(vl.x2, vl.y2) - - r1 := inRange(t1) - r2 := inRange(t2) - - switch { - case r1 && r2: - lin := linear(t[i], t[j]) - out = append(out, lin(t1), lin(t2)) - return out - - case !r1 && !r2: // the hole edge - out = append(out, t[i], t[j]) - return out - case !r1: - if t1 < 0 { - // below first -> clip by first - lin := linear(t[i], t[j]) - out = append(out, t[i], lin(t2)) - } else { - // above second -> clip by second - lin := linear(t[i], t[j]) - out = append(out, lin(t2), t[j]) - } - return out - case !r2: - if t2 < 0 { - // below first -> clip by first - lin := linear(t[i], t[j]) - out = append(out, t[i], lin(t1)) - } else { - // above second -> clip by second - lin := linear(t[i], t[j]) - out = append(out, lin(t1), t[j]) - } - return out - } - - case o1: - t1 := relative(t[i], t[j])(vl.x1, vl.y1) - if !inRange(t1) { - continue edges - } - out = append(out, linear(t[i], t[j])(t1)) - - case o2: - t2 := relative(t[i], t[j])(vl.x2, vl.y2) - if !inRange(t2) { - continue edges - } - out = append(out, linear(t[i], t[j])(t2)) - - default: - x, y, intersects := vl.line.Intersection(edge) - if !intersects { - continue edges - } - - // log.Println("Intersection -----------------------------") - t1 := relative(t[i], t[j])(x, y) - // log.Printf("relative pos: %f\n", t1) - if !inRange(t1) { - continue edges - } - - // log.Println("valid ***************") - - z := interpolate(t[j].Z, t[i].Z)(t1) - n := Vertex{x, y, z} - - if math.Signbit(s1) != math.Signbit(s2) { - // log.Println("\ton different sides") - // different sides -> vertex on edge. - out = append(out, n) - } else { // same side -> inside. Find peer - if len(out) > 0 { // we already have our peer - out = append(out, n) - return out - } - - for k := 0; k < 3; k++ { - if k == i { - continue - } - l := (k + 1) % 3 - other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y) - xo, yo, intersects := vl.line.Intersection(other) - if !intersects { - continue - } - t2 := relative(t[k], t[l])(xo, yo) - if !inRange(t2) { - continue - } - zo := interpolate(t[k].Z, t[l].Z)(t2) - - m := Vertex{xo, yo, zo} - - var xn, yn, xf, yf float64 - - // Which is nearer to current edge? - if math.Abs(s1) < math.Abs(s2) { - xn, yn = vl.x1, vl.y1 - xf, yf = vl.x2, vl.y2 - } else { - xn, yn = vl.x2, vl.y2 - xf, yf = vl.x1, vl.y1 - } - - if onPlane(other.Eval(xn, yn)) { - // triangle intersect in only point - // treat as no intersection - break edges - } - - pos := relative(n, m) - - tzn := pos(xn, yn) - tzf := pos(xf, yf) - - if !inRange(tzn) { - // if near is out of range far is, too. - return out - } - - lin := interpolate(n.Z, m.Z) - - if inRange(tzf) { - m.X = xf - m.Y = yf - m.Z = lin(tzf) - } // else m is clipping - - n.X = xn - n.Y = yn - n.Z = lin(tzn) - - out = append(out, n, m) - return out - } - } - } - } - - // supress single point touches. - if len(out) == 1 { - out = out[:0] - } - - return out -} - -// AsWKB returns a WKB representation of the given point cloud. -func (mpz MultiPointZ) AsWKB() []byte { - size := 1 + 4 + 4 + len(mpz)*(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(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 { - 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)) - } - - return buf.Bytes() -} - -func (mpz *MultiPointZ) 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.MultiPointZ: - return fmt.Errorf("unknown geometry type %x", geomType) - } - - var numPoints uint32 - if err := binary.Read(r, order, &numPoints); err != nil { - return err - } - - points := make(MultiPointZ, numPoints) - - for i := range points { - endian, err = r.ReadByte() - - 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) - } - - err = binary.Read(r, order, &geomType) - - switch { - case err != nil: - return err - case geomType != wkb.PointZ: - return fmt.Errorf("unknown geometry type %x", geomType) - } - - var x, y, z uint64 - if err = binary.Read(r, order, &x); err != nil { - return err - } - if err = binary.Read(r, order, &y); err != nil { - return err - } - if err = binary.Read(r, order, &z); err != nil { - return err - } - points[i] = Vertex{ - X: math.Float64frombits(x), - Y: math.Float64frombits(y), - Z: math.Float64frombits(z), - } - } - - *mpz = points - - return nil -}