Mercurial > gemma
view cmd/srsimplify/main.go @ 3770:71164b817d6e simplify-sounding-results
Calculate distance from planes.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 25 Jun 2019 21:04:16 +0200 |
parents | 6838526df94c |
children | b7530ed07561 |
line wrap: on
line source
// 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 main import ( "bufio" "flag" "fmt" "io" "log" "math" "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 simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ { if len(points) < 2 { return points } 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 - 2*tolerance xMin := min.X - 1 xMax := max.X + 1 yMin := min.Y - 1 yMax := max.Y + 1 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)) var rest int 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 } } rest++ } for i, part := range parts { log.Printf("%d: %d %.5f\n", i, len(part), maxDists[i]) } log.Printf("rest: %d\n", rest) // TODO: Implement me! return points } 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 = simplify(points, tolerance) if err := storeXYZ(points, os.Stdout); err != nil { log.Fatalf("err: %v\n", err) } }