Mercurial > gemma
changeset 4557:17cba8b447a6 iso-areas
Throw away own broken tracing algorithm.
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Tue, 01 Oct 2019 17:39:58 +0200 |
| parents | 04eba9dc917d |
| children | 07f632cd2625 |
| files | cmd/isoareas/algo.go cmd/isoareas/main.go |
| diffstat | 2 files changed, 0 insertions(+), 528 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/isoareas/algo.go Tue Oct 01 17:37:54 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,419 +0,0 @@ -// This is Free Software under GNU Affero General Public License v >= 3.0 -// without warranty, see README.md and license for details. -// -// SPDX-License-Identifier: AGPL-3.0-or-later -// License-Filename: LICENSES/AGPL-3.0.txt -// -// Copyright (C) 2019 by via donau -// – Österreichische Wasserstraßen-Gesellschaft mbH -// Software engineering by Intevation GmbH -// -// Author(s): -// * Sascha L. Teichmann <sascha.teichmann@intevation.de> - -package main - -import ( - "bufio" - "container/list" - "log" - "math" - "math/bits" - "os" - "runtime" - "sync" - - "gemma.intevation.de/gemma/pkg/octree" -) - -type indexedArc struct { - arc int32 - index int32 -} - -func glue(a, b octree.LineStringZ) octree.LineStringZ { - - a1, a2 := a[0], a[len(a)-1] - b1, b2 := b[0], b[len(b)-1] - - 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.EpsEquals2D(b1) { - c := make(octree.LineStringZ, len(a)-1+len(b)) - copy(c, a) - copy(c[len(a):], b[1:]) - return c - } - - if a1.EpsEquals2D(b1) { - c := make(octree.LineStringZ, len(a)-1+len(b)) - copy(c, b) - c[:len(b)].Reverse() - copy(c[len(b):], a[1:]) - return c - } - - if a2.EpsEquals2D(b2) { - c := make(octree.LineStringZ, len(a)-1+len(b)) - copy(c, a) - copy(c[len(a):], b[:len(b)-1]) - c[len(a):].Reverse() - return c - } - - return nil -} - -func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) { - - unique := map[int32]struct{}{} - for i := range cuts { - unique[cuts[i].arc] = struct{}{} - } - before := len(unique) - - origLen := int32(len(*arcs)) - - for i := range cuts { - if cuts[i].arc >= origLen { - // already has a connected arc assigned. - continue - } - - line := (*arcs)[cuts[i].arc] - - var modified []int - visited := map[int32]struct{}{} - - var recursive func(int32) - recursive = func(idx int32) { - visited[idx] = struct{}{} - - ns := neighbors(tri, idx) - for _, n := range ns { - n /= 3 - if _, already := visited[n]; already { - continue - } - - arcIdx := findArc(n, cuts) - if arcIdx == -1 { - continue - } - arc := cuts[arcIdx].arc - if arc >= origLen { - // already a new arc. - continue - } - oline := (*arcs)[arc] - - nline := glue(line, oline) - if nline == nil { - // not joinable - continue - } - line = nline - modified = append(modified, arcIdx) - recursive(cuts[arcIdx].index) - } - } - recursive(cuts[i].index) - if len(modified) > 0 { - // alloc a new line an fix the references. - nidx := int32(len(*arcs)) - *arcs = append(*arcs, line) - cuts[i].arc = nidx - for _, j := range modified { - cuts[j].arc = nidx - } - } - } - - unique = map[int32]struct{}{} - for i := range cuts { - unique[cuts[i].arc] = struct{}{} - } - log.Printf("unique arcs: before: %d after %d\n", - before, len(unique)) -} - -func intersectTriangles( - tri *octree.Triangulation, - heights []float64, - min, max octree.Vertex, -) [][]octree.LineStringZ { - - cuts := make([][]indexedArc, len(heights)) - - classes := make([][]int32, len(heights)+1) - - var arcs []octree.LineStringZ - -all: - for i := int32(0); i < int32(len(tri.Triangles)); i += 3 { - ti := tri.Triangles[i : i+3] - v0 := tri.Points[ti[0]] - v1 := tri.Points[ti[1]] - v2 := tri.Points[ti[2]] - min := math.Min(v0.Z, math.Min(v1.Z, v2.Z)) - - if min > heights[len(heights)-1] { - classes[len(classes)-1] = append(classes[len(classes)-1], i/3) - continue - } - max := math.Max(v0.Z, math.Max(v1.Z, v2.Z)) - - for j, h := range heights { - - var l float64 - if j > 0 { - l = heights[j-1] - } else { - l = -math.MaxFloat64 - } - - if l < min && max < h { - classes[j] = append(classes[j], i/3) - continue all - } - if min > h || max < h { - continue - } - t := octree.Triangle{v0, v1, v2} - cut := t.IntersectHorizontal(h) - if len(cut) >= 2 { - arc := int32(len(arcs)) - arcs = append(arcs, cut) - cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3}) - } - } - } - - // connect the arcs in a cut list to longer arcs. - - for _, c := range cuts { - 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)) - - jobs := make(chan int) - - var wg sync.WaitGroup - - worker := func() { - defer wg.Done() - - pb := polygonBuilder{open: list.New()} - - for i := range jobs { - usedArcs := map[int32]struct{}{} - var dupes int - var isolated, inside, found int - c := classes[i] - - 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 - - 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 - } - - 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 - - // 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) - } - } // 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) - /* - 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() - } - } - - 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)) - } - - return result - - // TODO: sew segments together. - -} - -type polygonBuilder struct { - polygons []octree.LineStringZ - - open *list.List -} - -func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) { - polygon := octree.LineStringZ{v0, v1, v2, v0} - 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] -} - -func findArc(needle int32, haystack []indexedArc) int { - lo, hi := 0, len(haystack)-1 - for lo <= hi { - mid := (hi-lo)/2 + lo - switch v := haystack[mid].index; { - case v < needle: - lo = mid + 1 - case v > needle: - hi = mid - 1 - default: - return mid - } - } - return -1 -} - -func contains(needle int32, haystack []int32) bool { - lo, hi := 0, len(haystack)-1 - for lo <= hi { - mid := (hi-lo)/2 + lo - switch v := haystack[mid]; { - case v < needle: - lo = mid + 1 - case v > needle: - hi = mid - 1 - default: - return true - } - } - return false -} - -func where(neighbors, indices []int32) byte { - var mask byte - for i, n := range neighbors { - if n >= 0 && contains(n/3, indices) { - mask |= 1 << i - } - } - return mask -}
--- a/cmd/isoareas/main.go Tue Oct 01 17:37:54 2019 +0200 +++ b/cmd/isoareas/main.go Tue Oct 01 17:39:58 2019 +0200 @@ -20,7 +20,6 @@ "io" "log" "math" - "math/rand" "os" "runtime" "strconv" @@ -169,114 +168,6 @@ return } -func dumpPolygonsSVG( - w io.Writer, - min, max octree.Vertex, - classPolygons [][]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 _, polygons := range classPolygons { - - r := byte(rnd.Intn(256)) - g := byte(rnd.Intn(256)) - b := byte(rnd.Intn(256)) - style := fmt.Sprintf("fill:#%02x%02x%02x", r, g, b) - - for _, polygon := range polygons { - - x := make([]int, len(polygon)) - y := make([]int, len(polygon)) - for i, p := range polygon { - x[i] = int(math.Round(px(p.X))) - y[i] = int(math.Round(py(p.Y))) - } - - canvas.Polygon(x, y, style) - } - - } - - canvas.End() - - 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() { cellSize := float64(0.5)