Mercurial > gemma
view cmd/isoareas/algo.go @ 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 | |
children |
line wrap: on
line source
// 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 }