Mercurial > gemma

comparison cmd/isoareas/main.go @ 4550:aa2d0006e742 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Write iso lines between classes to SVG for debugging.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 30 Sep 2019 16:34:55 +0200
parents 9c65cef72753
children b5b23b6d76f1
comparison
equal deleted inserted replaced
4549:9c65cef72753 4550:aa2d0006e742
21 "log" 21 "log"
22 "math" 22 "math"
23 "math/bits" 23 "math/bits"
24 "math/rand" 24 "math/rand"
25 "os" 25 "os"
26 "runtime"
26 "strconv" 27 "strconv"
27 "strings" 28 "strings"
29 "sync"
28 "time" 30 "time"
29 31
30 svg "github.com/ajstarks/svgo" 32 svg "github.com/ajstarks/svgo"
31 33
32 "gemma.intevation.de/gemma/pkg/octree" 34 "gemma.intevation.de/gemma/pkg/octree"
97 func glue(a, b octree.LineStringZ) octree.LineStringZ { 99 func glue(a, b octree.LineStringZ) octree.LineStringZ {
98 100
99 a1, a2 := a[0], a[len(a)-1] 101 a1, a2 := a[0], a[len(a)-1]
100 b1, b2 := b[0], b[len(b)-1] 102 b1, b2 := b[0], b[len(b)-1]
101 103
102 if a1.EpsEquals(b2) { 104 if a1.EpsEquals2D(b2) {
103 c := make(octree.LineStringZ, len(a)-1+len(b)) 105 c := make(octree.LineStringZ, len(a)-1+len(b))
104 copy(c, b) 106 copy(c, b)
105 copy(c[len(b):], a[1:]) 107 copy(c[len(b):], a[1:])
106 return c 108 return c
107 } 109 }
108 110
109 if a2.EpsEquals(b1) { 111 if a2.EpsEquals2D(b1) {
110 c := make(octree.LineStringZ, len(a)-1+len(b)) 112 c := make(octree.LineStringZ, len(a)-1+len(b))
111 copy(c, a) 113 copy(c, a)
112 copy(c[len(a):], b[1:]) 114 copy(c[len(a):], b[1:])
113 return c 115 return c
114 } 116 }
115 117
116 if a1.EpsEquals(b1) { 118 if a1.EpsEquals2D(b1) {
117 c := make(octree.LineStringZ, len(a)-1+len(b)) 119 c := make(octree.LineStringZ, len(a)-1+len(b))
118 copy(c, b) 120 copy(c, b)
119 c[:len(b)].Reverse() 121 c[:len(b)].Reverse()
120 copy(c[len(b):], a[1:]) 122 copy(c[len(b):], a[1:])
121 return c 123 return c
122 } 124 }
123 125
124 if a2.EpsEquals(b2) { 126 if a2.EpsEquals2D(b2) {
125 c := make(octree.LineStringZ, len(a)-1+len(b)) 127 c := make(octree.LineStringZ, len(a)-1+len(b))
126 copy(c, a) 128 copy(c, a)
127 copy(c[len(a):], b[:len(b)-1]) 129 copy(c[len(a):], b[:len(b)-1])
128 c[len(a):].Reverse() 130 c[len(a):].Reverse()
129 return c 131 return c
203 } 205 }
204 log.Printf("unique arcs: before: %d after %d\n", 206 log.Printf("unique arcs: before: %d after %d\n",
205 before, len(unique)) 207 before, len(unique))
206 } 208 }
207 209
208 func intersectTriangles(tri *octree.Triangulation, heights []float64) [][]octree.LineStringZ { 210 func intersectTriangles(
211 tri *octree.Triangulation,
212 heights []float64,
213 min, max octree.Vertex,
214 ) [][]octree.LineStringZ {
209 215
210 cuts := make([][]indexedArc, len(heights)) 216 cuts := make([][]indexedArc, len(heights))
211 217
212 classes := make([][]int32, len(heights)+1) 218 classes := make([][]int32, len(heights)+1)
213 219
257 263
258 for _, c := range cuts { 264 for _, c := range cuts {
259 connectArcs(tri, c, &arcs) 265 connectArcs(tri, c, &arcs)
260 } 266 }
261 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
262 result := make([][]octree.LineStringZ, len(classes)) 289 result := make([][]octree.LineStringZ, len(classes))
263 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
264 log.Println("inside classes:") 404 log.Println("inside classes:")
265 for i, c := range classes { 405 for i := range classes {
266 406 jobs <- i
267 pb := polygonBuilder{open: list.New()} 407 }
268 408 close(jobs)
269 usedArcs := map[int32]struct{}{} 409
270 var dupes int 410 wg.Wait()
271
272 var isolated, inside, found int
273 allInClass:
274 for _, idx := range c {
275 ns := neighbors(tri, idx)
276 mask := where(ns, c)
277 switch bits.OnesCount8(mask) {
278 case 3:
279 // Totally insides do not contribute to the geometries.
280 inside++
281 continue allInClass
282 case 0:
283 // Isolated are areas w/o connections to the rest.
284 isolated++
285 ti := tri.Triangles[idx*3 : idx*3+3]
286 pb.addTriangle(
287 tri.Points[ti[0]],
288 tri.Points[ti[1]],
289 tri.Points[ti[2]])
290 continue allInClass
291 default:
292 ti := tri.Triangles[idx*3 : idx*3+3]
293 for j := 0; j < 3; j++ {
294 if (mask & (1 << j)) == 0 {
295
296 var curr octree.LineStringZ
297
298 for l := i - 1; l <= i; l++ {
299 if l < 0 || l >= len(cuts) {
300 continue
301 }
302 arcIdx := findArc(ns[j]/3, cuts[l])
303 if arcIdx == -1 {
304 continue
305 }
306 found++
307 aIdx := cuts[l][arcIdx].arc
308 if _, already := usedArcs[aIdx]; already {
309 dupes++
310 continue
311 }
312 usedArcs[aIdx] = struct{}{}
313
314 curr = arcs[aIdx]
315 break
316 }
317
318 if curr == nil {
319 k := (j + 1) % 3
320 front := tri.Points[ti[j]]
321 back := tri.Points[ti[k]]
322 curr = octree.LineStringZ{front, back}
323 }
324
325 segment:
326 for e := pb.open.Front(); e != nil; e = e.Next() {
327 line := e.Value.(octree.LineStringZ)
328 if nline := glue(curr, line); nline != nil {
329 curr = nline
330 pb.open.Remove(e)
331 goto segment
332 }
333 } // all open
334
335 // check if closed
336 if curr[0].EpsEquals(curr[len(curr)-1]) {
337 if !curr.CCW() {
338 curr.Reverse()
339 }
340 pb.polygons = append(pb.polygons, curr)
341 } else {
342 pb.open.PushFront(curr)
343 }
344 } // for all border parts
345 }
346 }
347 }
348
349 log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n",
350 i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)
351
352 /*
353 for e := pb.open.Front(); e != nil; e = e.Next() {
354 line := e.Value.(octree.LineStringZ)
355 if !line.CCW() {
356 line.Reverse()
357 }
358 pb.polygons = append(pb.polygons, line)
359 }
360 */
361
362 result[i] = pb.polygons
363 }
364 411
365 log.Println("cuts:") 412 log.Println("cuts:")
366 for i, c := range cuts { 413 for i, c := range cuts {
367 log.Printf("\t%.3f: %d\n", heights[i], len(c)) 414 log.Printf("\t%.3f: %d\n", heights[i], len(c))
368 } 415 }
380 } 427 }
381 428
382 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) { 429 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
383 polygon := octree.LineStringZ{v0, v1, v2, v0} 430 polygon := octree.LineStringZ{v0, v1, v2, v0}
384 pb.polygons = append(pb.polygons, polygon) 431 pb.polygons = append(pb.polygons, polygon)
432 }
433
434 func (pb *polygonBuilder) reset() {
435 pb.polygons = pb.polygons[:0]
436 pb.open.Init()
385 } 437 }
386 438
387 func neighbors(t *octree.Triangulation, idx int32) []int32 { 439 func neighbors(t *octree.Triangulation, idx int32) []int32 {
388 idx *= 3 440 idx *= 3
389 return t.Halfedges[idx : idx+3] 441 return t.Halfedges[idx : idx+3]
500 canvas.End() 552 canvas.End()
501 553
502 return nil 554 return nil
503 } 555 }
504 556
557 func dumpArcsSVG(
558 w io.Writer,
559 min, max octree.Vertex,
560 cuts [][]indexedArc,
561 arcs []octree.LineStringZ,
562 ) (err error) {
563
564 ratio := (max.X - min.X) / (max.Y - min.Y)
565
566 log.Printf("ratio: %.2f\n", ratio)
567
568 const width = 50000
569 height := int(math.Ceil(width * ratio))
570
571 px := linear(min.X, 0, max.X, width)
572 py := linear(min.Y, float64(height), max.Y, 0)
573
574 out := bufio.NewWriter(w)
575 defer func() { err = out.Flush() }()
576
577 canvas := svg.New(out)
578 canvas.Start(width, height)
579
580 rnd := rand.New(rand.NewSource(42))
581
582 for _, cut := range cuts {
583
584 usedArcs := map[int32]struct{}{}
585
586 r := byte(rnd.Intn(256))
587 g := byte(rnd.Intn(256))
588 b := byte(rnd.Intn(256))
589
590 style := fmt.Sprintf("fill:none;stroke:#%02x%02x%02x", r, g, b)
591
592 for i := range cut {
593 idx := cut[i].arc
594 if _, already := usedArcs[idx]; already {
595 continue
596 }
597 usedArcs[idx] = struct{}{}
598 arc := arcs[idx]
599 x := make([]int, len(arc))
600 y := make([]int, len(arc))
601 for i, p := range arc {
602 x[i] = int(math.Round(px(p.X)))
603 y[i] = int(math.Round(py(p.Y)))
604 }
605
606 canvas.Polyline(x, y, style)
607 }
608 }
609
610 canvas.End()
611
612 return nil
613 }
614
505 func main() { 615 func main() {
506 616
507 heights, err := octree.ParseClassBreaks(classBreaks) 617 heights, err := octree.ParseClassBreaks(classBreaks)
508 check(err) 618 check(err)
509 log.Printf("num classes: %d\n", len(heights)) 619 log.Printf("num classes: %d\n", len(heights))
527 log.Printf("triangulation took %v\n", time.Since(start)) 637 log.Printf("triangulation took %v\n", time.Since(start))
528 638
529 log.Printf("number of triangles: %d\n", len(tri.Triangles)/3) 639 log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)
530 640
531 start = time.Now() 641 start = time.Now()
532 polygons := intersectTriangles(tri, heights) 642 polygons := intersectTriangles(tri, heights, min, max)
533 log.Printf("intersecting triangles took %v.\n", time.Since(start)) 643 log.Printf("intersecting triangles took %v.\n", time.Since(start))
534 644
535 check(dumpPolygonsSVG(os.Stdout, min, max, polygons)) 645 check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
536 } 646 }