Mercurial > gemma
view cmd/isoareas/main.go @ 4536:3130c005abef iso-areas
Measure time of triangulation correctly.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 25 Sep 2019 17:28:13 +0200 |
parents | 508075a5694e |
children | 46c9a731a686 |
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" "sort" "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 int32 } cuts := make([][]indexedCut, len(heights)) classes := make([][]int32, len(heights)+1) all: for i := int32(0); i < int32(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)) if min > heights[len(heights)-1] { classes[len(classes)-1] = append(classes[len(classes)-1], i/3) continue } max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) for j, h := range heights { var l float64 if j > 0 { l = heights[j-1] } else { l = -math.MaxFloat64 } if l < min && max < h { classes[j] = append(classes[j], i/3) continue all } if min > h || max < h { continue } cut := t.IntersectHorizontal(h) if len(cut) >= 2 { cuts[j] = append(cuts[j], indexedCut{cut, i / 3}) } } } log.Println("inside classes:") for i, c := range classes { remain := make([]int32, 0, len(c)) for _, idx := range c { // Eliminate triangles that do not contribute // any border to the final outline. if !inner(tri, idx, c) { remain = append(remain, idx) } } log.Printf("\t%d (%d) %.2f%%\n", len(remain), len(c), 100*float64(len(remain))/float64(len(c))) classes[i] = remain } log.Println("cuts:") for i, c := range cuts { log.Printf("\t%.3f: %d\n", heights[i], len(c)) } // TODO: sew segments together. } func neighbors(t *octree.Triangulation, idx int32) []int32 { idx *= 3 return t.Halfedges[idx : idx+3] } func inner(t *octree.Triangulation, idx int32, indices []int32) bool { for _, n := range neighbors(t, idx) { if n < 0 { return false } n /= 3 if p := sort.Search(len(indices), func(i int) bool { return indices[i] >= n }); p >= len(indices) || indices[p] != n { return false } } return true } func main() { heights, err := octree.ParseClassBreaks(classBreaks) check(err) log.Printf("num classes: %d\n", len(heights)) start := time.Now() xyz, err := loadXYZ(os.Stdin) check(err) log.Printf("loading took %v\n", time.Since(start)) 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) start = time.Now() tri, err := octree.Triangulate(xyz) check(err) log.Printf("triangulation took %v\n", time.Since(start)) log.Printf("number of triangles: %d\n", len(tri.Triangles)/3) start = time.Now() intersectTriangles(tri, heights) log.Printf("intersecting triangles took %v.\n", time.Since(start)) }