Mercurial > gemma

comparison cmd/isoareas/main.go @ 4533:3998a9ab69c6 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Find out which triangles are fully inside classes.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 25 Sep 2019 13:04:08 +0200
parents 769f723c2581
children b7d31a449dd2
comparison
equal deleted inserted replaced
4532:769f723c2581 4533:3998a9ab69c6
85 85
86 func intersectTriangles(tri *octree.Triangulation, heights []float64) { 86 func intersectTriangles(tri *octree.Triangulation, heights []float64) {
87 87
88 type indexedCut struct { 88 type indexedCut struct {
89 cut octree.LineStringZ 89 cut octree.LineStringZ
90 index int 90 index int32
91 } 91 }
92 92
93 cuts := make([][]indexedCut, len(heights)) 93 cuts := make([][]indexedCut, len(heights))
94 94
95 for i := 0; i < len(tri.Triangles); i += 3 { 95 classes := make([][]int32, len(heights)+1)
96
97 all:
98 for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
96 ti := tri.Triangles[i : i+3] 99 ti := tri.Triangles[i : i+3]
97 t := octree.Triangle{ 100 t := octree.Triangle{
98 tri.Points[ti[0]], 101 tri.Points[ti[0]],
99 tri.Points[ti[1]], 102 tri.Points[ti[1]],
100 tri.Points[ti[2]], 103 tri.Points[ti[2]],
101 } 104 }
102 min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) 105 min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z))
106
107 if min > heights[len(heights)-1] {
108 classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
109 continue
110 }
103 max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) 111 max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z))
104 112
105 for j, h := range heights { 113 for j, h := range heights {
106 if h < min { 114
115 var l float64
116 if j > 0 {
117 l = heights[j-1]
118 } else {
119 l = -math.MaxFloat64
120 }
121
122 if l < min && max < h {
123 classes[j] = append(classes[j], i/3)
124 continue all
125 }
126 if min > h || max < h {
107 continue 127 continue
108 } 128 }
109 if h > max { 129 cut := t.IntersectHorizontal(h)
110 break 130 if len(cut) >= 2 {
131 cuts[j] = append(cuts[j], indexedCut{cut, i / 3})
111 } 132 }
112 cut := t.IntersectHorizontal(h)
113 if len(cut) < 2 {
114 continue
115 }
116 cuts[j] = append(cuts[j], indexedCut{cut, i})
117 } 133 }
118 } 134 }
119 135
120 log.Println("cuts") 136 log.Println("inside classes:")
121 for i := range cuts { 137 for _, c := range classes {
122 log.Printf("%.3f: %d\n", heights[i], len(cuts[i])) 138 log.Printf("\t%d\n", len(c))
139 }
140
141 log.Println("cuts:")
142 for i, c := range cuts {
143 log.Printf("\t%.3f: %d\n", heights[i], len(c))
123 } 144 }
124 145
125 // TODO: sew segments together. 146 // TODO: sew segments together.
126 147
127 } 148 }
128 149
129 func main() { 150 func main() {
130 151
131 heights, err := octree.ParseClassBreaks(classBreaks) 152 heights, err := octree.ParseClassBreaks(classBreaks)
132 check(err) 153 check(err)
154 log.Printf("num classes: %d\n", len(heights))
133 155
156 start := time.Now()
134 xyz, err := loadXYZ(os.Stdin) 157 xyz, err := loadXYZ(os.Stdin)
135 check(err) 158 check(err)
159 log.Printf("loading took %v\n", time.Since(start))
136 160
137 log.Printf("num vertices: %d\n", len(xyz)) 161 log.Printf("num vertices: %d\n", len(xyz))
138 162
139 min, max := minMax(xyz) 163 min, max := minMax(xyz)
140 164
141 log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n", 165 log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
142 min.X, min.Y, min.Z, 166 min.X, min.Y, min.Z,
143 max.X, max.Y, max.Z) 167 max.X, max.Y, max.Z)
144 168
145 heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z) 169 start = time.Now()
146
147 log.Printf("classes: %d\n", len(heights))
148
149 tri, err := octree.Triangulate(xyz) 170 tri, err := octree.Triangulate(xyz)
150 check(err) 171 check(err)
172 log.Printf("triangulation took %v\n", time.Since(time.Now()))
151 173
152 start := time.Now() 174 log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)
175
176 start = time.Now()
153 intersectTriangles(tri, heights) 177 intersectTriangles(tri, heights)
154 log.Printf("intersecting triangles took %v.\n", time.Since(start)) 178 log.Printf("intersecting triangles took %v.\n", time.Since(start))
155 } 179 }