Mercurial > gemma

changeset 4550:aa2d0006e742 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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).