# HG changeset patch # User Sascha L. Teichmann # Date 1561972247 -7200 # Node ID bb62c98fcf054900663f02973e187189ae7758d2 # Parent fd6d62b08af78b2838c29db1477a9b3b484a29c6# Parent 57c8cfff4920776ae68510c2ba7633e3e275e3a4 Merged simplify-sounding-results into default. diff -r fd6d62b08af7 -r bb62c98fcf05 cmd/srsimplify/main.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/srsimplify/main.go Mon Jul 01 11:10:47 2019 +0200 @@ -0,0 +1,91 @@ +// 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 + +package main + +import ( + "bufio" + "flag" + "fmt" + "io" + "log" + "os" + "strconv" + "strings" + + "gemma.intevation.de/gemma/pkg/octree" +) + +func loadXYZ(r io.Reader) (octree.MultiPointZ, error) { + + scanner := bufio.NewScanner(r) + + points := make(octree.MultiPointZ, 0, 2000000) + + var x, y, z float64 + var err error + + for scanner.Scan() { + line := strings.TrimSpace(scanner.Text()) + if len(line) == 0 || strings.HasPrefix(line, "#") { + continue + } + parts := strings.SplitN(line, " ", 3) + if len(parts) != 3 { + continue + } + + if x, err = strconv.ParseFloat(parts[0], 64); err != nil { + return nil, err + } + if y, err = strconv.ParseFloat(parts[1], 64); err != nil { + return nil, err + } + if z, err = strconv.ParseFloat(parts[2], 64); err != nil { + return nil, err + } + points = append(points, octree.Vertex{X: x, Y: y, Z: z}) + } + + return points, nil +} + +func storeXYZ(points octree.MultiPointZ, w io.Writer) error { + out := bufio.NewWriter(w) + for i := range points { + fmt.Fprintf(out, "%.5f,%.5f,%.5f\n", + points[i].X, points[i].Y, points[i].Z) + } + return out.Flush() +} + +func main() { + + var tolerance float64 + + flag.Float64Var(&tolerance, "t", 0.1, "accepted tolerance (shorthand)") + flag.Float64Var(&tolerance, "tolerance", 0.1, "accepted tolerance") + + flag.Parse() + + points, err := loadXYZ(os.Stdin) + if err != nil { + log.Fatalf("err: %v\n", err) + } + + points = points.Simplify(tolerance) + + if err := storeXYZ(points, os.Stdout); err != nil { + log.Fatalf("err: %v\n", err) + } +} diff -r fd6d62b08af7 -r bb62c98fcf05 pkg/octree/simplify.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/simplify.go Mon Jul 01 11:10:47 2019 +0200 @@ -0,0 +1,189 @@ +// 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 + +package octree + +import ( + "math" +) + +func (points MultiPointZ) Simplify(tolerance float64) MultiPointZ { + + if len(points) < 2 { + return points + } + + if tolerance < 0 { + tolerance = -tolerance + } + + min := Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64} + max := Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64} + + var maxIdx int + + for i, v := range points { + min.Minimize(v) + + if v.X < min.X { + min.X = v.X + } + if v.X > max.X { + max.X = v.X + } + if v.Y < min.Y { + min.Y = v.Y + } + if v.Y > max.Y { + max.Y = v.Y + } + if v.Z < min.Z { + min.Z = v.Z + } + if v.Z > max.Z { + max.Z = v.Z + maxIdx = i + } + + max.Maximize(v) + } + + /* + log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n", + min.X, min.Y, min.Z, + max.X, max.Y, max.Z) + */ + + below := min.Z - 3*tolerance + xMin := min.X - tolerance + xMax := max.X + tolerance + yMin := min.Y - tolerance + yMax := max.Y + tolerance + + corners := []Vertex{ + {xMin, yMin, below}, + {xMax, yMin, below}, + {xMax, yMax, below}, + {xMin, yMax, below}, + } + + top := points[maxIdx] + + tris := make([]Triangle, len(corners)) + planes := make([]Plane3D, len(corners)) + + for i, v1 := range corners { + v2 := corners[(i+1)%len(corners)] + tris[i] = Triangle{v1, v2, top} + planes[i] = tris[i].Plane3D() + } + + parts := make([][]Vertex, len(tris)) + + maxDists := make([]float64, len(planes)) + maxIdxs := make([]int, len(planes)) + +nextPoint: + for i, v := range points { + if i == maxIdx { + continue + } + + for j := range tris { + if tris[j].Contains(v.X, v.Y) { + if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] { + maxDists[j] = dist + maxIdxs[j] = len(parts[j]) + } + parts[j] = append(parts[j], v) + continue nextPoint + } + } + } + + result := make([]Vertex, 0, len(points)) + + var handleTriangle func(*Triangle, float64, int, []Vertex) bool + + handleTriangle = func( + t *Triangle, + maxDist float64, maxIdx int, + points []Vertex, + ) bool { + if maxDist <= tolerance { + return false + } + + if len(points) == 1 { + result = append(result, points[0]) + } + + var ( + tris [3]Triangle + planes [3]Plane3D + maxDists [3]float64 + maxIdxs [3]int + parts [3][]Vertex + ) + + top := points[maxIdx] + for i := 0; i < 3; i++ { + tris[i] = Triangle{t[i], t[(i+1)%3], top} + planes[i] = tris[i].Plane3D() + } + + nextPoint: + for i, v := range points { + if i == maxIdx { + continue + } + + for j := range tris { + if tris[j].Contains(v.X, v.Y) { + if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] { + maxDists[j] = dist + maxIdxs[j] = len(parts[j]) + } + parts[j] = append(parts[j], v) + continue nextPoint + } + } + } + + var found bool + for i, part := range parts { + if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) { + found = true + } + } + + if found { + result = append(result, top) + } + + return found + } + + var found bool + for i, part := range parts { + if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) { + found = true + } + } + + if found { + result = append(result, top) + } + + return result +} diff -r fd6d62b08af7 -r bb62c98fcf05 pkg/octree/vertex.go --- a/pkg/octree/vertex.go Thu Jun 27 16:27:48 2019 +0200 +++ b/pkg/octree/vertex.go Mon Jul 01 11:10:47 2019 +0200 @@ -110,6 +110,10 @@ 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{