Mercurial > gemma

comparison cmd/isoareas/main.go @ 4551:b5b23b6d76f1 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move own algorith to separate file.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 01 Oct 2019 11:07:33 +0200
parents aa2d0006e742
children 7ed5a4a94efc
comparison
equal deleted inserted replaced
4550:aa2d0006e742 4551:b5b23b6d76f1
13 13
14 package main 14 package main
15 15
16 import ( 16 import (
17 "bufio" 17 "bufio"
18 "container/list"
19 "fmt" 18 "fmt"
20 "io" 19 "io"
21 "log" 20 "log"
22 "math" 21 "math"
23 "math/bits"
24 "math/rand" 22 "math/rand"
25 "os" 23 "os"
26 "runtime"
27 "strconv" 24 "strconv"
28 "strings" 25 "strings"
29 "sync"
30 "time" 26 "time"
31 27
32 svg "github.com/ajstarks/svgo" 28 svg "github.com/ajstarks/svgo"
33 29
34 "gemma.intevation.de/gemma/pkg/octree" 30 "gemma.intevation.de/gemma/pkg/octree"
87 83
88 func check(err error) { 84 func check(err error) {
89 if err != nil { 85 if err != nil {
90 log.Fatalf("error: %v\n", err) 86 log.Fatalf("error: %v\n", err)
91 } 87 }
92 }
93
94 type indexedArc struct {
95 arc int32
96 index int32
97 }
98
99 func glue(a, b octree.LineStringZ) octree.LineStringZ {
100
101 a1, a2 := a[0], a[len(a)-1]
102 b1, b2 := b[0], b[len(b)-1]
103
104 if a1.EpsEquals2D(b2) {
105 c := make(octree.LineStringZ, len(a)-1+len(b))
106 copy(c, b)
107 copy(c[len(b):], a[1:])
108 return c
109 }
110
111 if a2.EpsEquals2D(b1) {
112 c := make(octree.LineStringZ, len(a)-1+len(b))
113 copy(c, a)
114 copy(c[len(a):], b[1:])
115 return c
116 }
117
118 if a1.EpsEquals2D(b1) {
119 c := make(octree.LineStringZ, len(a)-1+len(b))
120 copy(c, b)
121 c[:len(b)].Reverse()
122 copy(c[len(b):], a[1:])
123 return c
124 }
125
126 if a2.EpsEquals2D(b2) {
127 c := make(octree.LineStringZ, len(a)-1+len(b))
128 copy(c, a)
129 copy(c[len(a):], b[:len(b)-1])
130 c[len(a):].Reverse()
131 return c
132 }
133
134 return nil
135 }
136
137 func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) {
138
139 unique := map[int32]struct{}{}
140 for i := range cuts {
141 unique[cuts[i].arc] = struct{}{}
142 }
143 before := len(unique)
144
145 origLen := int32(len(*arcs))
146
147 for i := range cuts {
148 if cuts[i].arc >= origLen {
149 // already has a connected arc assigned.
150 continue
151 }
152
153 line := (*arcs)[cuts[i].arc]
154
155 var modified []int
156 visited := map[int32]struct{}{}
157
158 var recursive func(int32)
159 recursive = func(idx int32) {
160 visited[idx] = struct{}{}
161
162 ns := neighbors(tri, idx)
163 for _, n := range ns {
164 n /= 3
165 if _, already := visited[n]; already {
166 continue
167 }
168
169 arcIdx := findArc(n, cuts)
170 if arcIdx == -1 {
171 continue
172 }
173 arc := cuts[arcIdx].arc
174 if arc >= origLen {
175 // already a new arc.
176 continue
177 }
178 oline := (*arcs)[arc]
179
180 nline := glue(line, oline)
181 if nline == nil {
182 // not joinable
183 continue
184 }
185 line = nline
186 modified = append(modified, arcIdx)
187 recursive(cuts[arcIdx].index)
188 }
189 }
190 recursive(cuts[i].index)
191 if len(modified) > 0 {
192 // alloc a new line an fix the references.
193 nidx := int32(len(*arcs))
194 *arcs = append(*arcs, line)
195 cuts[i].arc = nidx
196 for _, j := range modified {
197 cuts[j].arc = nidx
198 }
199 }
200 }
201
202 unique = map[int32]struct{}{}
203 for i := range cuts {
204 unique[cuts[i].arc] = struct{}{}
205 }
206 log.Printf("unique arcs: before: %d after %d\n",
207 before, len(unique))
208 }
209
210 func intersectTriangles(
211 tri *octree.Triangulation,
212 heights []float64,
213 min, max octree.Vertex,
214 ) [][]octree.LineStringZ {
215
216 cuts := make([][]indexedArc, len(heights))
217
218 classes := make([][]int32, len(heights)+1)
219
220 var arcs []octree.LineStringZ
221
222 all:
223 for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
224 ti := tri.Triangles[i : i+3]
225 v0 := tri.Points[ti[0]]
226 v1 := tri.Points[ti[1]]
227 v2 := tri.Points[ti[2]]
228 min := math.Min(v0.Z, math.Min(v1.Z, v2.Z))
229
230 if min > heights[len(heights)-1] {
231 classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
232 continue
233 }
234 max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))
235
236 for j, h := range heights {
237
238 var l float64
239 if j > 0 {
240 l = heights[j-1]
241 } else {
242 l = -math.MaxFloat64
243 }
244
245 if l < min && max < h {
246 classes[j] = append(classes[j], i/3)
247 continue all
248 }
249 if min > h || max < h {
250 continue
251 }
252 t := octree.Triangle{v0, v1, v2}
253 cut := t.IntersectHorizontal(h)
254 if len(cut) >= 2 {
255 arc := int32(len(arcs))
256 arcs = append(arcs, cut)
257 cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3})
258 }
259 }
260 }
261
262 // connect the arcs in a cut list to longer arcs.
263
264 for _, c := range cuts {
265 connectArcs(tri, c, &arcs)
266 }
267
268 func() {
269 out, err := os.Create("arcs.svg")
270 if err != nil {
271 log.Printf("err: %v\n", err)
272 return
273 }
274 defer func() {
275 out.Close()
276 log.Println("writing arcs done")
277 }()
278
279 buf := bufio.NewWriter(out)
280
281 dumpArcsSVG(
282 buf,
283 min, max,
284 cuts,
285 arcs)
286 buf.Flush()
287 }()
288
289 result := make([][]octree.LineStringZ, len(classes))
290
291 jobs := make(chan int)
292
293 var wg sync.WaitGroup
294
295 worker := func() {
296 defer wg.Done()
297
298 pb := polygonBuilder{open: list.New()}
299
300 for i := range jobs {
301 usedArcs := map[int32]struct{}{}
302 var dupes int
303 var isolated, inside, found int
304 c := classes[i]
305
306 allInClass:
307 for _, idx := range c {
308 ns := neighbors(tri, idx)
309 mask := where(ns, c)
310 switch bits.OnesCount8(mask) {
311 case 3:
312 // Totally insides do not contribute to the geometries.
313 inside++
314 continue allInClass
315 case 0:
316 // Isolated are areas w/o connections to the rest.
317 isolated++
318 ti := tri.Triangles[idx*3 : idx*3+3]
319 pb.addTriangle(
320 tri.Points[ti[0]],
321 tri.Points[ti[1]],
322 tri.Points[ti[2]])
323 continue allInClass
324 default:
325 ti := tri.Triangles[idx*3 : idx*3+3]
326 for j := 0; j < 3; j++ {
327 if (mask & (1 << j)) == 0 {
328
329 var curr octree.LineStringZ
330
331 for l := i - 1; l <= i; l++ {
332 if l < 0 || l >= len(cuts) {
333 continue
334 }
335 arcIdx := findArc(ns[j]/3, cuts[l])
336 if arcIdx == -1 {
337 continue
338 }
339 found++
340 aIdx := cuts[l][arcIdx].arc
341 if _, already := usedArcs[aIdx]; already {
342 dupes++
343 continue
344 }
345 usedArcs[aIdx] = struct{}{}
346
347 curr = arcs[aIdx]
348 break
349 }
350
351 if curr == nil {
352 k := (j + 1) % 3
353 front := tri.Points[ti[j]]
354 back := tri.Points[ti[k]]
355 curr = octree.LineStringZ{front, back}
356 }
357
358 segment:
359 for e := pb.open.Front(); e != nil; e = e.Next() {
360 line := e.Value.(octree.LineStringZ)
361 if nline := glue(curr, line); nline != nil {
362 curr = nline
363 pb.open.Remove(e)
364 goto segment
365 }
366 } // all open
367
368 // check if closed
369 if curr[0].EpsEquals2D(curr[len(curr)-1]) {
370 if !curr.CCW() {
371 curr.Reverse()
372 }
373 pb.polygons = append(pb.polygons, curr)
374 } else {
375 pb.open.PushFront(curr)
376 }
377 } // for all border parts
378 }
379 }
380 }
381
382 log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n",
383 i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)
384 /*
385 for e := pb.open.Front(); e != nil; e = e.Next() {
386 line := e.Value.(octree.LineStringZ)
387 if !line.CCW() {
388 line.Reverse()
389 }
390 pb.polygons = append(pb.polygons, line)
391 }
392 */
393
394 result[i] = pb.polygons
395 pb.reset()
396 }
397 }
398
399 for n := runtime.NumCPU(); n >= 0; n-- {
400 wg.Add(1)
401 go worker()
402 }
403
404 log.Println("inside classes:")
405 for i := range classes {
406 jobs <- i
407 }
408 close(jobs)
409
410 wg.Wait()
411
412 log.Println("cuts:")
413 for i, c := range cuts {
414 log.Printf("\t%.3f: %d\n", heights[i], len(c))
415 }
416
417 return result
418
419 // TODO: sew segments together.
420
421 }
422
423 type polygonBuilder struct {
424 polygons []octree.LineStringZ
425
426 open *list.List
427 }
428
429 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
430 polygon := octree.LineStringZ{v0, v1, v2, v0}
431 pb.polygons = append(pb.polygons, polygon)
432 }
433
434 func (pb *polygonBuilder) reset() {
435 pb.polygons = pb.polygons[:0]
436 pb.open.Init()
437 }
438
439 func neighbors(t *octree.Triangulation, idx int32) []int32 {
440 idx *= 3
441 return t.Halfedges[idx : idx+3]
442 }
443
444 func findArc(needle int32, haystack []indexedArc) int {
445 lo, hi := 0, len(haystack)-1
446 for lo <= hi {
447 mid := (hi-lo)/2 + lo
448 switch v := haystack[mid].index; {
449 case v < needle:
450 lo = mid + 1
451 case v > needle:
452 hi = mid - 1
453 default:
454 return mid
455 }
456 }
457 return -1
458 }
459
460 func contains(needle int32, haystack []int32) bool {
461 lo, hi := 0, len(haystack)-1
462 for lo <= hi {
463 mid := (hi-lo)/2 + lo
464 switch v := haystack[mid]; {
465 case v < needle:
466 lo = mid + 1
467 case v > needle:
468 hi = mid - 1
469 default:
470 return true
471 }
472 }
473 return false
474 }
475
476 func where(neighbors, indices []int32) byte {
477 var mask byte
478 for i, n := range neighbors {
479 if n >= 0 && contains(n/3, indices) {
480 mask |= 1 << i
481 }
482 }
483 return mask
484 } 88 }
485 89
486 func linear(x1, y1, x2, y2 float64) func(float64) float64 { 90 func linear(x1, y1, x2, y2 float64) func(float64) float64 {
487 // f(x1) = y1 91 // f(x1) = y1
488 // f(x2) = y2 92 // f(x2) = y2