Mercurial > gemma
view cmd/isoareas/main.go @ 4544:e5ae16f6d846 iso-areas
Lay foundation to join cut segments to longer arcs.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 27 Sep 2019 13:21:20 +0200 |
parents | 2a707732331f |
children | e7970d84cb4f |
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" "fmt" "io" "log" "math" "math/bits" "math/rand" "os" "strconv" "strings" "time" svg "github.com/ajstarks/svgo" "gemma.intevation.de/gemma/pkg/octree" ) const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff` func loadXYZ(r io.Reader) (octree.MultiPointZ, error) { scanner := bufio.NewScanner(r) points := make(octree.MultiPointZ, 0, 2000000) var x, y, z float64 var err error for scanner.Scan() { line := strings.TrimSpace(scanner.Text()) if len(line) == 0 || strings.HasPrefix(line, "#") { continue } parts := strings.SplitN(line, " ", 3) if len(parts) != 3 { continue } if x, err = strconv.ParseFloat(parts[0], 64); err != nil { return nil, err } if y, err = strconv.ParseFloat(parts[1], 64); err != nil { return nil, err } if z, err = strconv.ParseFloat(parts[2], 64); err != nil { return nil, err } points = append(points, octree.Vertex{X: x, Y: y, Z: z}) } return points, nil } func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) { if len(points) == 0 { return octree.Vertex{}, octree.Vertex{} } min, max := points[0], points[0] for _, v := range points { min.Minimize(v) max.Maximize(v) } return min, max } func check(err error) { if err != nil { log.Fatalf("error: %v\n", err) } } type indexedArc struct { arc int32 index int32 } func glue(a, b octree.LineStringZ) (octree.LineStringZ, bool) { a1, a2 := a[0], a[len(a)-1] b1, b2 := b[0], b[len(b)-1] if a1.EpsEquals(b2) { c := make(octree.LineStringZ, len(a)-1+len(b)) copy(c, b) copy(c[len(b):], a[1:]) return c, true } if a2.EpsEquals(b1) { c := make(octree.LineStringZ, len(a)-1+len(b)) copy(c, a) copy(c[len(a):], b[1:]) return c, true } if a1.EpsEquals(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, true } if a2.EpsEquals(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, true } return nil, false } func connectArcs(tri *octree.Triangulation, cut []indexedArc, arcs *[]octree.LineStringZ) { // TODO: Implement me! } func intersectTriangles(tri *octree.Triangulation, heights []float64) [][]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, i / 3}) } } } // connect the arcs in a cut list to longer arcs. for _, c := range cuts { connectArcs(tri, c, &arcs) } result := make([][]octree.LineStringZ, len(classes)) log.Println("inside classes:") for i, c := range classes { pb := polygonBuilder{open: list.New()} var isolated, inside int allInClass: for _, idx := range c { ns := neighbors(tri, idx) mask := where(ns, c) switch bits.OnesCount8(mask) { case 0: // Totally insides do not contribute to the geometries. inside++ continue allInClass case 3: // 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 { // TODO: Look into cuts to see // if there are real intersections var curr octree.LineStringZ for l := i - 1; l <= i; l++ { if l < 0 || l >= len(cuts) { continue } if arc := findArc(ns[j], cuts[l]); arc >= 0 { curr = arcs[arc] 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) first, last := line[0], line[len(line)-1] front, back := curr[0], curr[len(curr)-1] if back.EpsEquals(first) { curr = append(curr, line[1:]...) pb.open.Remove(e) goto segment } if back.EpsEquals(last) { line.Reverse() curr = append(curr, line[1:]...) pb.open.Remove(e) goto segment } if front.EpsEquals(last) { curr = append(line, curr[1:]...) pb.open.Remove(e) goto segment } if front.EpsEquals(first) { line.Reverse() curr = append(line, curr[1:]...) pb.open.Remove(e) goto segment } } // all open // check if closed if curr[0].EpsEquals(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\n", i, inside, isolated, pb.open.Len(), len(pb.polygons)) result[i] = pb.polygons } 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 neighbors(t *octree.Triangulation, idx int32) []int32 { idx *= 3 return t.Halfedges[idx : idx+3] } func findArc(needle int32, haystack []indexedArc) int32 { 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 haystack[mid].arc } } 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 // y1 = x1*a + b <=> b = y1 - x1*a // y2 = x2*a + b // y1 - y2 = a*(x1 - x2) // a = (y1-y2)/(x1 - x2) for x1 != x2 if x1 == x2 { return func(float64) float64 { return 0.5 * (y1 + y2) } } a := (y1 - y2) / (x1 - x2) b := y1 - x1*a return func(x float64) float64 { return x*a + b } } 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 = 1000 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 main() { heights, err := octree.ParseClassBreaks(classBreaks) check(err) log.Printf("num classes: %d\n", len(heights)) start := time.Now() xyz, err := loadXYZ(os.Stdin) check(err) log.Printf("loading took %v\n", time.Since(start)) log.Printf("num vertices: %d\n", len(xyz)) min, max := minMax(xyz) log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n", min.X, min.Y, min.Z, max.X, max.Y, max.Z) start = time.Now() tri, err := octree.Triangulate(xyz) check(err) log.Printf("triangulation took %v\n", time.Since(start)) log.Printf("number of triangles: %d\n", len(tri.Triangles)/3) start = time.Now() polygons := intersectTriangles(tri, heights) log.Printf("intersecting triangles took %v.\n", time.Since(start)) check(dumpPolygonsSVG(os.Stdout, min, max, polygons)) }