Mercurial > gemma

annotate cmd/isoareas/main.go @ 4534:b7d31a449dd2 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Find out how many inner triangles can be skipped because they are totally inside a class polygon.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 25 Sep 2019 15:39:55 +0200
parents 3998a9ab69c6
children 508075a5694e
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"
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
22 "sort"
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
23 "strconv"
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
24 "strings"
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
25 "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
26
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
27 "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
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
30 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
31
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
32 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
33
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
34 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
35
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
36 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
37
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 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
39 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
40
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
41 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
42 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
43 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
44 continue
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
45 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
46 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
47 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
48 continue
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
51 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
52 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
53 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
54 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
55 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
56 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
57 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
58 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
59 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
60 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
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
63 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
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
66 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
67 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
68 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
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
71 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
72
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
73 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
74 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
75 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
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
78 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
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
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
81 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
82 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
83 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
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 }
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
86
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
87 func intersectTriangles(tri *octree.Triangulation, heights []float64) {
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
88
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
89 type indexedCut struct {
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
90 cut octree.LineStringZ
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
91 index int32
4532
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
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
94 cuts := make([][]indexedCut, len(heights))
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
95
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
96 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
97
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
98 all:
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
99 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
100 ti := tri.Triangles[i : i+3]
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
101 t := octree.Triangle{
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
102 tri.Points[ti[0]],
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
103 tri.Points[ti[1]],
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
104 tri.Points[ti[2]],
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
105 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
106 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
107
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
108 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
109 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
110 continue
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
111 }
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
112 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
113
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
114 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
115
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
116 var l float64
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
117 if j > 0 {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
118 l = heights[j-1]
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
119 } else {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
120 l = -math.MaxFloat64
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
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
123 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
124 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
125 continue all
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
126 }
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
127 if min > h || max < h {
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
128 continue
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
129 }
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
130 cut := t.IntersectHorizontal(h)
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
131 if len(cut) >= 2 {
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
132 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
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 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
136
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
137 log.Println("inside classes:")
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
138 for _, c := range classes {
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
139 var eliminated int
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 {
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
141 if inner(tri, idx, c) {
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
142 eliminated++
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
143 }
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
144 }
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 log.Printf("\t%d (%d) %.2f%%\n",
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 len(c)-eliminated, len(c),
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 100*float64(eliminated)/float64(len(c)))
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
148 }
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
149
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
150 log.Println("cuts:")
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
151 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
152 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
153 }
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
154
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
155 // TODO: sew segments together.
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
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
159 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
160 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
161 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
162 }
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
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 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
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 for _, n := range neighbors(t, idx) {
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
167 if n < 0 {
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
168 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
169 }
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
170 n /= 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
171 p := sort.Search(len(indices), func(i int) 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
172 return indices[i] >= n
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
173 })
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
174 if p >= len(indices) || indices[p] != n {
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
175 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
176 }
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
177 }
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
178
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
179 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
180 }
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
181
4531
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
182 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
183
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
184 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
185 check(err)
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
186 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
187
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
188 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
189 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
190 check(err)
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
191 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
192
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
193 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
194
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
195 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
196
c9b6be8d81c8 Started a branch to calculate the iso areas for a given triangulation.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents:
diff changeset
197 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
198 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
199 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
200
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
201 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
202 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
203 check(err)
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
204 log.Printf("triangulation took %v\n", time.Since(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
205
4533
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
206 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
207
3998a9ab69c6 Find out which triangles are fully inside classes.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4532
diff changeset
208 start = time.Now()
4532
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
209 intersectTriangles(tri, heights)
769f723c2581 Cut triangles against class breaks.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 4531
diff changeset
210 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
211 }