Mercurial > gemma
changeset 4551:b5b23b6d76f1 iso-areas
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 |
files | cmd/isoareas/algo.go cmd/isoareas/main.go |
diffstat | 2 files changed, 419 insertions(+), 396 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/isoareas/algo.go Tue Oct 01 11:07:33 2019 +0200 @@ -0,0 +1,419 @@ +// 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 Mon Sep 30 16:34:55 2019 +0200 +++ b/cmd/isoareas/main.go Tue Oct 01 11:07:33 2019 +0200 @@ -15,18 +15,14 @@ import ( "bufio" - "container/list" "fmt" "io" "log" "math" - "math/bits" "math/rand" "os" - "runtime" "strconv" "strings" - "sync" "time" svg "github.com/ajstarks/svgo" @@ -91,398 +87,6 @@ } } -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 -} - func linear(x1, y1, x2, y2 float64) func(float64) float64 { // f(x1) = y1 // f(x2) = y2