# HG changeset patch # User Sascha L. Teichmann # Date 1561972099 -7200 # Node ID 545304d3ff936f57ee2b9734fd7e24e2e411efa3 # Parent b7530ed075613d4cf0c6f796fb7b091e3f8f60ae Moved simplification algorithm to octree package. diff -r b7530ed07561 -r 545304d3ff93 cmd/srsimplify/main.go --- a/cmd/srsimplify/main.go Wed Jun 26 12:06:28 2019 +0200 +++ b/cmd/srsimplify/main.go Mon Jul 01 11:08:19 2019 +0200 @@ -19,7 +19,6 @@ "fmt" "io" "log" - "math" "os" "strconv" "strings" @@ -70,175 +69,6 @@ return out.Flush() } -func simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ { - - if len(points) < 2 { - return points - } - - if tolerance < 0 { - tolerance = -tolerance - } - - min := octree.Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64} - max := octree.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 := []octree.Vertex{ - {xMin, yMin, below}, - {xMax, yMin, below}, - {xMax, yMax, below}, - {xMin, yMax, below}, - } - - top := points[maxIdx] - - tris := make([]octree.Triangle, len(corners)) - planes := make([]octree.Plane3D, len(corners)) - - for i, v1 := range corners { - v2 := corners[(i+1)%len(corners)] - tris[i] = octree.Triangle{v1, v2, top} - planes[i] = tris[i].Plane3D() - } - - parts := make([][]octree.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([]octree.Vertex, 0, len(points)) - - var handleTriangle func(*octree.Triangle, float64, int, []octree.Vertex) bool - - handleTriangle = func( - t *octree.Triangle, - maxDist float64, maxIdx int, - points []octree.Vertex, - ) bool { - if maxDist <= tolerance { - return false - } - - if len(points) == 1 { - result = append(result, points[0]) - } - - var ( - tris [3]octree.Triangle - planes [3]octree.Plane3D - maxDists [3]float64 - maxIdxs [3]int - parts [3][]octree.Vertex - ) - - top := points[maxIdx] - for i := 0; i < 3; i++ { - tris[i] = octree.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 -} - func main() { var tolerance float64 @@ -253,7 +83,7 @@ log.Fatalf("err: %v\n", err) } - points = simplify(points, tolerance) + points = points.Simplify(tolerance) if err := storeXYZ(points, os.Stdout); err != nil { log.Fatalf("err: %v\n", err) diff -r b7530ed07561 -r 545304d3ff93 pkg/octree/simplify.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/simplify.go Mon Jul 01 11:08:19 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 +}