Mercurial > gemma
changeset 4550:aa2d0006e742 iso-areas
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 |
files | cmd/isoareas/main.go pkg/octree/vertex.go |
diffstat | 2 files changed, 207 insertions(+), 89 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/isoareas/main.go Mon Sep 30 14:59:20 2019 +0200 +++ b/cmd/isoareas/main.go Mon Sep 30 16:34:55 2019 +0200 @@ -23,8 +23,10 @@ "math/bits" "math/rand" "os" + "runtime" "strconv" "strings" + "sync" "time" svg "github.com/ajstarks/svgo" @@ -99,21 +101,21 @@ a1, a2 := a[0], a[len(a)-1] b1, b2 := b[0], b[len(b)-1] - if a1.EpsEquals(b2) { + if a1.EpsEquals2D(b2) { c := make(octree.LineStringZ, len(a)-1+len(b)) copy(c, b) copy(c[len(b):], a[1:]) return c } - if a2.EpsEquals(b1) { + if a2.EpsEquals2D(b1) { c := make(octree.LineStringZ, len(a)-1+len(b)) copy(c, a) copy(c[len(a):], b[1:]) return c } - if a1.EpsEquals(b1) { + if a1.EpsEquals2D(b1) { c := make(octree.LineStringZ, len(a)-1+len(b)) copy(c, b) c[:len(b)].Reverse() @@ -121,7 +123,7 @@ return c } - if a2.EpsEquals(b2) { + if a2.EpsEquals2D(b2) { c := make(octree.LineStringZ, len(a)-1+len(b)) copy(c, a) copy(c[len(a):], b[:len(b)-1]) @@ -205,7 +207,11 @@ before, len(unique)) } -func intersectTriangles(tri *octree.Triangulation, heights []float64) [][]octree.LineStringZ { +func intersectTriangles( + tri *octree.Triangulation, + heights []float64, + min, max octree.Vertex, +) [][]octree.LineStringZ { cuts := make([][]indexedArc, len(heights)) @@ -259,109 +265,150 @@ connectArcs(tri, c, &arcs) } + func() { + out, err := os.Create("arcs.svg") + if err != nil { + log.Printf("err: %v\n", err) + return + } + defer func() { + out.Close() + log.Println("writing arcs done") + }() + + buf := bufio.NewWriter(out) + + dumpArcsSVG( + buf, + min, max, + cuts, + arcs) + buf.Flush() + }() + result := make([][]octree.LineStringZ, len(classes)) - log.Println("inside classes:") - for i, c := range classes { + jobs := make(chan int) + + var wg sync.WaitGroup + + worker := func() { + defer wg.Done() pb := polygonBuilder{open: list.New()} - usedArcs := map[int32]struct{}{} - var dupes int + for i := range jobs { + usedArcs := map[int32]struct{}{} + var dupes int + var isolated, inside, found int + c := classes[i] - var isolated, inside, found int - allInClass: - for _, idx := range c { - ns := neighbors(tri, idx) - mask := where(ns, c) - switch bits.OnesCount8(mask) { - case 3: - // Totally insides do not contribute to the geometries. - inside++ - continue allInClass - case 0: - // Isolated are areas w/o connections to the rest. - isolated++ - ti := tri.Triangles[idx*3 : idx*3+3] - pb.addTriangle( - tri.Points[ti[0]], - tri.Points[ti[1]], - tri.Points[ti[2]]) - continue allInClass - default: - ti := tri.Triangles[idx*3 : idx*3+3] - for j := 0; j < 3; j++ { - if (mask & (1 << j)) == 0 { + allInClass: + for _, idx := range c { + ns := neighbors(tri, idx) + mask := where(ns, c) + switch bits.OnesCount8(mask) { + case 3: + // Totally insides do not contribute to the geometries. + inside++ + continue allInClass + case 0: + // Isolated are areas w/o connections to the rest. + isolated++ + ti := tri.Triangles[idx*3 : idx*3+3] + pb.addTriangle( + tri.Points[ti[0]], + tri.Points[ti[1]], + tri.Points[ti[2]]) + continue allInClass + default: + ti := tri.Triangles[idx*3 : idx*3+3] + for j := 0; j < 3; j++ { + if (mask & (1 << j)) == 0 { - var curr octree.LineStringZ + var curr octree.LineStringZ - for l := i - 1; l <= i; l++ { - if l < 0 || l >= len(cuts) { - continue - } - arcIdx := findArc(ns[j]/3, cuts[l]) - if arcIdx == -1 { - continue + for l := i - 1; l <= i; l++ { + if l < 0 || l >= len(cuts) { + continue + } + arcIdx := findArc(ns[j]/3, cuts[l]) + if arcIdx == -1 { + continue + } + found++ + aIdx := cuts[l][arcIdx].arc + if _, already := usedArcs[aIdx]; already { + dupes++ + continue + } + usedArcs[aIdx] = struct{}{} + + curr = arcs[aIdx] + break } - found++ - aIdx := cuts[l][arcIdx].arc - if _, already := usedArcs[aIdx]; already { - dupes++ - continue + + if curr == nil { + k := (j + 1) % 3 + front := tri.Points[ti[j]] + back := tri.Points[ti[k]] + curr = octree.LineStringZ{front, back} } - usedArcs[aIdx] = struct{}{} - - curr = arcs[aIdx] - break - } - - if curr == nil { - k := (j + 1) % 3 - front := tri.Points[ti[j]] - back := tri.Points[ti[k]] - curr = octree.LineStringZ{front, back} - } - segment: - for e := pb.open.Front(); e != nil; e = e.Next() { - line := e.Value.(octree.LineStringZ) - if nline := glue(curr, line); nline != nil { - curr = nline - pb.open.Remove(e) - goto segment - } - } // all open + segment: + for e := pb.open.Front(); e != nil; e = e.Next() { + line := e.Value.(octree.LineStringZ) + if nline := glue(curr, line); nline != nil { + curr = nline + pb.open.Remove(e) + goto segment + } + } // all open - // check if closed - if curr[0].EpsEquals(curr[len(curr)-1]) { - if !curr.CCW() { - curr.Reverse() + // check if closed + if curr[0].EpsEquals2D(curr[len(curr)-1]) { + if !curr.CCW() { + curr.Reverse() + } + pb.polygons = append(pb.polygons, curr) + } else { + pb.open.PushFront(curr) } - pb.polygons = append(pb.polygons, curr) - } else { - pb.open.PushFront(curr) - } - } // for all border parts + } // for all border parts + } } } - } - log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n", - i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found) + log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n", + i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found) + /* + for e := pb.open.Front(); e != nil; e = e.Next() { + line := e.Value.(octree.LineStringZ) + if !line.CCW() { + line.Reverse() + } + pb.polygons = append(pb.polygons, line) + } + */ - /* - for e := pb.open.Front(); e != nil; e = e.Next() { - line := e.Value.(octree.LineStringZ) - if !line.CCW() { - line.Reverse() - } - pb.polygons = append(pb.polygons, line) - } - */ + result[i] = pb.polygons + pb.reset() + } + } - result[i] = pb.polygons + for n := runtime.NumCPU(); n >= 0; n-- { + wg.Add(1) + go worker() } + log.Println("inside classes:") + for i := range classes { + jobs <- i + } + close(jobs) + + wg.Wait() + log.Println("cuts:") for i, c := range cuts { log.Printf("\t%.3f: %d\n", heights[i], len(c)) @@ -384,6 +431,11 @@ pb.polygons = append(pb.polygons, polygon) } +func (pb *polygonBuilder) reset() { + pb.polygons = pb.polygons[:0] + pb.open.Init() +} + func neighbors(t *octree.Triangulation, idx int32) []int32 { idx *= 3 return t.Halfedges[idx : idx+3] @@ -502,6 +554,64 @@ return nil } +func dumpArcsSVG( + w io.Writer, + min, max octree.Vertex, + cuts [][]indexedArc, + arcs []octree.LineStringZ, +) (err error) { + + ratio := (max.X - min.X) / (max.Y - min.Y) + + log.Printf("ratio: %.2f\n", ratio) + + const width = 50000 + height := int(math.Ceil(width * ratio)) + + px := linear(min.X, 0, max.X, width) + py := linear(min.Y, float64(height), max.Y, 0) + + out := bufio.NewWriter(w) + defer func() { err = out.Flush() }() + + canvas := svg.New(out) + canvas.Start(width, height) + + rnd := rand.New(rand.NewSource(42)) + + for _, cut := range cuts { + + usedArcs := map[int32]struct{}{} + + r := byte(rnd.Intn(256)) + g := byte(rnd.Intn(256)) + b := byte(rnd.Intn(256)) + + style := fmt.Sprintf("fill:none;stroke:#%02x%02x%02x", r, g, b) + + for i := range cut { + idx := cut[i].arc + if _, already := usedArcs[idx]; already { + continue + } + usedArcs[idx] = struct{}{} + arc := arcs[idx] + x := make([]int, len(arc)) + y := make([]int, len(arc)) + for i, p := range arc { + x[i] = int(math.Round(px(p.X))) + y[i] = int(math.Round(py(p.Y))) + } + + canvas.Polyline(x, y, style) + } + } + + canvas.End() + + return nil +} + func main() { heights, err := octree.ParseClassBreaks(classBreaks) @@ -529,7 +639,7 @@ log.Printf("number of triangles: %d\n", len(tri.Triangles)/3) start = time.Now() - polygons := intersectTriangles(tri, heights) + polygons := intersectTriangles(tri, heights, min, max) log.Printf("intersecting triangles took %v.\n", time.Since(start)) check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
--- a/pkg/octree/vertex.go Mon Sep 30 14:59:20 2019 +0200 +++ b/pkg/octree/vertex.go Mon Sep 30 16:34:55 2019 +0200 @@ -483,6 +483,14 @@ math.Abs(v.Y-w.Y) < eps && math.Abs(v.Z-w.Z) < eps } +// EpsEquals returns true if v and w are equal component-wise +// in the X/Y plane with the values within a epsilon range. +func (v Vertex) EpsEquals2D(w Vertex) bool { + const eps = 1e-5 + return math.Abs(v.X-w.X) < eps && + math.Abs(v.Y-w.Y) < eps +} + // JoinOnLine joins the the elements of a given multi line string // under the assumption that the segments are all on the line segment // from (x1, y1) to (x2, y2).