# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1537435923 -7200
# Node ID be90ab542aa738c42cc617f562195be4d9424d65
# Parent  a2f107f1e4e71351639270951a5a9f395541835e
octree: contouring: Do the math to calculate the intersection points of the triangles and the planes.

diff -r a2f107f1e4e7 -r be90ab542aa7 cmd/octree2contour/loader.go
--- a/cmd/octree2contour/loader.go	Thu Sep 20 00:13:47 2018 +0200
+++ b/cmd/octree2contour/loader.go	Thu Sep 20 11:32:03 2018 +0200
@@ -129,7 +129,7 @@
 	return tree, nil
 }
 
-func (ot *octree) horizontal(h float64, fn func([]int32)) {
+func (ot *octree) horizontal(h float64, fn func(*triangle)) {
 
 	type frame struct {
 		pos int32
@@ -178,14 +178,16 @@
 					continue
 				}
 				tri := ot.triangles[idx]
-				v0 := ot.vertices[tri[0]]
-				v1 := ot.vertices[tri[1]]
-				v2 := ot.vertices[tri[2]]
+				t := triangle{
+					ot.vertices[tri[0]],
+					ot.vertices[tri[1]],
+					ot.vertices[tri[2]],
+				}
 
-				if !(math.Min(v0.z, math.Min(v1.z, v2.z)) > h ||
-					math.Max(v0.z, math.Max(v1.z, v2.z)) < h) {
+				if !(math.Min(t[0].z, math.Min(t[1].z, t[2].z)) > h ||
+					math.Max(t[0].z, math.Max(t[1].z, t[2].z)) < h) {
 					dupes[idx] = struct{}{}
-					fn(tri)
+					fn(&t)
 				}
 			}
 		}
diff -r a2f107f1e4e7 -r be90ab542aa7 cmd/octree2contour/main.go
--- a/cmd/octree2contour/main.go	Thu Sep 20 00:13:47 2018 +0200
+++ b/cmd/octree2contour/main.go	Thu Sep 20 11:32:03 2018 +0200
@@ -17,7 +17,7 @@
 
 type result struct {
 	h         float64
-	triangles int
+	numPoints int
 }
 
 func processLevels(
@@ -28,21 +28,23 @@
 ) {
 	defer wg.Done()
 	for h := range jobs {
-		var triangles int
-		tree.horizontal(h, func([]int32) {
-			triangles++
+		var points []vertex
+		tree.horizontal(h, func(t *triangle) {
+			points = t.intersectH(h, points)
 		})
-		results <- result{h, triangles}
+		results <- result{h, len(points)}
 	}
 }
 
 func process(tree *octree) {
-	var triangles int
+	var numPoints int
 	start := time.Now()
 	if *one {
-		tree.horizontal(*step, func([]int32) {
-			triangles++
+		var points []vertex
+		tree.horizontal(*step, func(t *triangle) {
+			points = t.intersectH(*step, points)
 		})
+		numPoints = len(points)
 	} else {
 
 		results := make(chan result)
@@ -78,14 +80,13 @@
 			return all[i].h < all[j].h
 		})
 
-		var sum int
 		for i := range all {
 			a := &all[i]
-			sum += a.triangles
-			log.Printf("level %f: %d\n", a.h, a.triangles)
+			numPoints += a.numPoints
+			log.Printf("level %f: %d\n", a.h, a.numPoints)
 		}
-		log.Printf("num triangles: %d\n", sum)
 	}
+	log.Printf("num points: %d\n", numPoints)
 	log.Printf("processing took: %s\n", time.Since(start))
 }
 
diff -r a2f107f1e4e7 -r be90ab542aa7 cmd/octree2contour/vertex.go
--- a/cmd/octree2contour/vertex.go	Thu Sep 20 00:13:47 2018 +0200
+++ b/cmd/octree2contour/vertex.go	Thu Sep 20 11:32:03 2018 +0200
@@ -38,6 +38,22 @@
 	}
 }
 
+func (v vertex) add(w vertex) vertex {
+	return vertex{
+		v.x + w.x,
+		v.y + w.y,
+		v.z + w.z,
+	}
+}
+
+func (v vertex) scale(s float64) vertex {
+	return vertex{
+		s * v.x,
+		s * v.y,
+		s * v.z,
+	}
+}
+
 func interpolate(v1, v2 vertex) func(vertex) vertex {
 	v2 = v2.sub(v1)
 	return func(s vertex) vertex {
@@ -52,3 +68,67 @@
 func (a vertex) less(b vertex) bool {
 	return a.x < b.x || a.y < b.y || a.z < b.z
 }
+
+type line [2]vertex
+
+func newLine(p1, p2 vertex) line {
+	return line{
+		p2.sub(p1),
+		p1,
+	}
+}
+
+func (l line) eval(t float64) vertex {
+	return l[0].scale(t).add(l[1])
+}
+
+func (l line) intersectH(h float64) vertex {
+	t := (h - l[1].z) / l[0].z
+	return l.eval(t)
+}
+
+type triangle [3]vertex
+
+func side(z, h float64) int {
+	switch {
+	case z < h:
+		return -1
+	case z > h:
+		return +1
+	}
+	return 0
+}
+
+func (t *triangle) intersectH(h float64, points []vertex) []vertex {
+	sides := [3]int{
+		side(t[0].z, h),
+		side(t[1].z, h),
+		side(t[2].z, h),
+	}
+
+	for i := 0; i < 3; i++ {
+		j := (i + 1) % 3
+		si := sides[i]
+		sj := sides[j]
+
+		switch {
+		case si == 0 && sj == 0:
+			// both on plane
+			points = append(points, t[i], t[j])
+		case si == 0 && sj != 0:
+			// first on plane
+			points = append(points, t[i])
+		case si != 0 && sj == 0:
+			// second on plane
+			points = append(points, t[j])
+		case si == sj:
+			// both on same side
+		default:
+			// real intersection
+			v := newLine(t[i], t[j]).intersectH(h)
+			points = append(points, v)
+		}
+	}
+
+	return points
+}