Mercurial > gemma

annotate cmd/isoareas/main.go @ 4539:99928231b679 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Code simplification.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 26 Sep 2019 10:59:17 +0200
parents 46c9a731a686
children e01b6e7cbe7d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
1 // This is Free Software under GNU Affero General Public License v >= 3.0
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
2 // without warranty, see README.md and license for details.
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
3 //
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
4 // SPDX-License-Identifier: AGPL-3.0-or-later
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
5 // License-Filename: LICENSES/AGPL-3.0.txt
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
6 //
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
7 // Copyright (C) 2019 by via donau
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
8 // – Österreichische Wasserstraßen-Gesellschaft mbH
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
9 // Software engineering by Intevation GmbH
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
10 //
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
11 // Author(s):
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
12 // * Sascha L. Teichmann <sascha.teichmann@intevation.de>
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
13
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
14 package main
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
15
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
16 import (
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
17 "bufio"
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
18 "io"
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
19 "log"
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
20 "math"
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
21 "os"
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
22 "strconv"
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
23 "strings"
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
24 "time"
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
25
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
26 "gemma.intevation.de/gemma/pkg/octree"
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
27 )
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
28
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
29 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`
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
30
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
31 func loadXYZ(r io.Reader) (octree.MultiPointZ, error) {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
32
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
33 scanner := bufio.NewScanner(r)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
34
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
35 points := make(octree.MultiPointZ, 0, 2000000)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
36
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
37 var x, y, z float64
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
38 var err error
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
39
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
40 for scanner.Scan() {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
41 line := strings.TrimSpace(scanner.Text())
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
42 if len(line) == 0 || strings.HasPrefix(line, "#") {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
43 continue
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
44 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
45 parts := strings.SplitN(line, " ", 3)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
46 if len(parts) != 3 {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
47 continue
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
48 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
49
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
50 if x, err = strconv.ParseFloat(parts[0], 64); err != nil {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
51 return nil, err
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
52 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
53 if y, err = strconv.ParseFloat(parts[1], 64); err != nil {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
54 return nil, err
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
55 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
56 if z, err = strconv.ParseFloat(parts[2], 64); err != nil {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
57 return nil, err
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
58 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
59 points = append(points, octree.Vertex{X: x, Y: y, Z: z})
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
60 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
61
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
62 return points, nil
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
63 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
64
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
65 func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
66 if len(points) == 0 {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
67 return octree.Vertex{}, octree.Vertex{}
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
68 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
69
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
70 min, max := points[0], points[0]
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
71
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
72 for _, v := range points {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
73 min.Minimize(v)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
74 max.Maximize(v)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
75 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
76
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
77 return min, max
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
78 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
79
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
80 func check(err error) {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
81 if err != nil {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
82 log.Fatalf("error: %v\n", err)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
83 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
84 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
85
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
86 func intersectTriangles(tri *octree.Triangulation, heights []float64) {
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
87
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
88 type indexedCut struct {
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
89 cut octree.LineStringZ
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
90 index int32
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
91 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
92
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
93 cuts := make([][]indexedCut, len(heights))
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
94
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
95 classes := make([][]int32, len(heights)+1)
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
96
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
97 all:
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
98 for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
99 ti := tri.Triangles[i : i+3]
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
100 t := octree.Triangle{
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
101 tri.Points[ti[0]],
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
102 tri.Points[ti[1]],
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
103 tri.Points[ti[2]],
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
104 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
105 min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z))
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
106
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
107 if min > heights[len(heights)-1] {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
108 classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
109 continue
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
110 }
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
111 max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z))
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
112
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
113 for j, h := range heights {
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
114
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
115 var l float64
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
116 if j > 0 {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
117 l = heights[j-1]
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
118 } else {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
119 l = -math.MaxFloat64
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
120 }
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
121
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
122 if l < min && max < h {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
123 classes[j] = append(classes[j], i/3)
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
124 continue all
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
125 }
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
126 if min > h || max < h {
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
127 continue
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
128 }
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
129 cut := t.IntersectHorizontal(h)
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
130 if len(cut) >= 2 {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
131 cuts[j] = append(cuts[j], indexedCut{cut, i / 3})
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
132 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
133 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
134 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
135
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
136 log.Println("inside classes:")
4535
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
137 for i, c := range classes {
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
138 remain := make([]int32, 0, len(c))
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
139
4534
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
140 for _, idx := range c {
4535
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
141 // Eliminate triangles that do not contribute
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
142 // any border to the final outline.
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
143 if !inner(tri, idx, c) {
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
144 remain = append(remain, idx)
4534
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
145 }
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
146 }
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
147 log.Printf("\t%d (%d) %.2f%%\n",
4535
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
148 len(remain), len(c),
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
149 100*float64(len(remain))/float64(len(c)))
508075a5694e Only keep the indices of the triangles of the border of the iso patches.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4534
diff changeset
150 classes[i] = remain
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
151 }
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
152
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
153 log.Println("cuts:")
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
154 for i, c := range cuts {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
155 log.Printf("\t%.3f: %d\n", heights[i], len(c))
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
156 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
157
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
158 // TODO: sew segments together.
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
159
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
160 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
161
4534
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
162 func neighbors(t *octree.Triangulation, idx int32) []int32 {
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
163 idx *= 3
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
164 return t.Halfedges[idx : idx+3]
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
165 }
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
166
4538
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
167 func contains(needle int32, haystack []int32) bool {
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
168 lo, hi := 0, len(haystack)-1
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
169 for lo <= hi {
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
170 mid := (hi-lo)/2 + lo
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
171 switch v := haystack[mid]; {
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
172 case v < needle:
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
173 lo = mid + 1
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
174 case v > needle:
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
175 hi = mid - 1
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
176 default:
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
177 return true
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
178 }
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
179 }
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
180 return false
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
181 }
46c9a731a686 Use custom binary search in triangle indices.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4536
diff changeset
182
4534
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
183 func inner(t *octree.Triangulation, idx int32, indices []int32) bool {
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
184
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
185 for _, n := range neighbors(t, idx) {
4539
99928231b679 Code simplification.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4538
diff changeset
186 if n < 0 || !contains(n/3, indices) {
4534
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
187 return false
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
188 }
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
189 }
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
190
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
191 return true
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
192 }
b7d31a449dd2 Find out how many inner triangles can be skipped because they are totally inside a class polygon.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4533
diff changeset
193
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
194 func main() {
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
195
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
196 heights, err := octree.ParseClassBreaks(classBreaks)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
197 check(err)
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
198 log.Printf("num classes: %d\n", len(heights))
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
199
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
200 start := time.Now()
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
201 xyz, err := loadXYZ(os.Stdin)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
202 check(err)
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
203 log.Printf("loading took %v\n", time.Since(start))
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
204
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
205 log.Printf("num vertices: %d\n", len(xyz))
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
206
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
207 min, max := minMax(xyz)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
208
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
209 log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
210 min.X, min.Y, min.Z,
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
211 max.X, max.Y, max.Z)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
212
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
213 start = time.Now()
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
214 tri, err := octree.Triangulate(xyz)
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
215 check(err)
4536
3130c005abef Measure time of triangulation correctly.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4535
diff changeset
216 log.Printf("triangulation took %v\n", time.Since(start))
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
217
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
218 log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
219
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
220 start = time.Now()
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
221 intersectTriangles(tri, heights)
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
222 log.Printf("intersecting triangles took %v.\n", time.Since(start))
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
223 }