Mercurial > gemma

comparison cmd/isoareas/main.go @ 4532:769f723c2581 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
4531:c9b6be8d81c8 4532:769f723c2581
15 15
16 import ( 16 import (
17 "bufio" 17 "bufio"
18 "io" 18 "io"
19 "log" 19 "log"
20 "math"
20 "os" 21 "os"
21 "strconv" 22 "strconv"
22 "strings" 23 "strings"
24 "time"
23 25
24 "gemma.intevation.de/gemma/pkg/octree" 26 "gemma.intevation.de/gemma/pkg/octree"
25 ) 27 )
26 28
27 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` 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`
79 if err != nil { 81 if err != nil {
80 log.Fatalf("error: %v\n", err) 82 log.Fatalf("error: %v\n", err)
81 } 83 }
82 } 84 }
83 85
86 func intersectTriangles(tri *octree.Triangulation, heights []float64) {
87
88 type indexedCut struct {
89 cut octree.LineStringZ
90 index int
91 }
92
93 cuts := make([][]indexedCut, len(heights))
94
95 for i := 0; i < len(tri.Triangles); i += 3 {
96 ti := tri.Triangles[i : i+3]
97 t := octree.Triangle{
98 tri.Points[ti[0]],
99 tri.Points[ti[1]],
100 tri.Points[ti[2]],
101 }
102 min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z))
103 max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z))
104
105 for j, h := range heights {
106 if h < min {
107 continue
108 }
109 if h > max {
110 break
111 }
112 cut := t.IntersectHorizontal(h)
113 if len(cut) < 2 {
114 continue
115 }
116 cuts[j] = append(cuts[j], indexedCut{cut, i})
117 }
118 }
119
120 log.Println("cuts")
121 for i := range cuts {
122 log.Printf("%.3f: %d\n", heights[i], len(cuts[i]))
123 }
124
125 // TODO: sew segments together.
126
127 }
128
84 func main() { 129 func main() {
85 130
86 heights, err := octree.ParseClassBreaks(classBreaks) 131 heights, err := octree.ParseClassBreaks(classBreaks)
87 check(err) 132 check(err)
88 133
97 min.X, min.Y, min.Z, 142 min.X, min.Y, min.Z,
98 max.X, max.Y, max.Z) 143 max.X, max.Y, max.Z)
99 144
100 heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z) 145 heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z)
101 146
102 log.Println("classes:") 147 log.Printf("classes: %d\n", len(heights))
103 for i, v := range heights {
104 log.Printf("\t%d: %.3f\n", i+1, v)
105 }
106 148
107 tri, err := octree.Triangulate(xyz) 149 tri, err := octree.Triangulate(xyz)
108 check(err) 150 check(err)
109 151
110 _ = tri 152 start := time.Now()
153 intersectTriangles(tri, heights)
154 log.Printf("intersecting triangles took %v.\n", time.Since(start))
111 } 155 }