Mercurial > gemma
view cmd/isoareas/main.go @ 4532:769f723c2581 iso-areas
Cut triangles against class breaks.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 24 Sep 2019 18:03:43 +0200 |
parents | c9b6be8d81c8 |
children | 3998a9ab69c6 |
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" "io" "log" "math" "os" "strconv" "strings" "time" "gemma.intevation.de/gemma/pkg/octree" ) const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff` 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 minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) { if len(points) == 0 { return octree.Vertex{}, octree.Vertex{} } min, max := points[0], points[0] for _, v := range points { min.Minimize(v) max.Maximize(v) } return min, max } func check(err error) { if err != nil { log.Fatalf("error: %v\n", err) } } func intersectTriangles(tri *octree.Triangulation, heights []float64) { type indexedCut struct { cut octree.LineStringZ index int } cuts := make([][]indexedCut, len(heights)) for i := 0; i < len(tri.Triangles); i += 3 { ti := tri.Triangles[i : i+3] t := octree.Triangle{ tri.Points[ti[0]], tri.Points[ti[1]], tri.Points[ti[2]], } min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) for j, h := range heights { if h < min { continue } if h > max { break } cut := t.IntersectHorizontal(h) if len(cut) < 2 { continue } cuts[j] = append(cuts[j], indexedCut{cut, i}) } } log.Println("cuts") for i := range cuts { log.Printf("%.3f: %d\n", heights[i], len(cuts[i])) } // TODO: sew segments together. } func main() { heights, err := octree.ParseClassBreaks(classBreaks) check(err) xyz, err := loadXYZ(os.Stdin) check(err) log.Printf("num vertices: %d\n", len(xyz)) min, max := minMax(xyz) log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n", min.X, min.Y, min.Z, max.X, max.Y, max.Z) heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z) log.Printf("classes: %d\n", len(heights)) tri, err := octree.Triangulate(xyz) check(err) start := time.Now() intersectTriangles(tri, heights) log.Printf("intersecting triangles took %v.\n", time.Since(start)) }