# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1569409448 -7200
# Node ID 3998a9ab69c65b98149c9f8eb8f768828574f4a3
# Parent  769f723c2581039b4daabc1f83cda289c2164118
Find out which triangles are fully inside classes.

diff -r 769f723c2581 -r 3998a9ab69c6 cmd/isoareas/main.go
--- a/cmd/isoareas/main.go	Tue Sep 24 18:03:43 2019 +0200
+++ b/cmd/isoareas/main.go	Wed Sep 25 13:04:08 2019 +0200
@@ -87,12 +87,15 @@
 
 	type indexedCut struct {
 		cut   octree.LineStringZ
-		index int
+		index int32
 	}
 
 	cuts := make([][]indexedCut, len(heights))
 
-	for i := 0; i < len(tri.Triangles); i += 3 {
+	classes := make([][]int32, len(heights)+1)
+
+all:
+	for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
 		ti := tri.Triangles[i : i+3]
 		t := octree.Triangle{
 			tri.Points[ti[0]],
@@ -100,26 +103,44 @@
 			tri.Points[ti[2]],
 		}
 		min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z))
+
+		if min > heights[len(heights)-1] {
+			classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
+			continue
+		}
 		max := math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z))
 
 		for j, h := range heights {
-			if h < min {
+
+			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
 			}
-			if h > max {
-				break
+			cut := t.IntersectHorizontal(h)
+			if len(cut) >= 2 {
+				cuts[j] = append(cuts[j], indexedCut{cut, i / 3})
 			}
-			cut := t.IntersectHorizontal(h)
-			if len(cut) < 2 {
-				continue
-			}
-			cuts[j] = append(cuts[j], indexedCut{cut, i})
 		}
 	}
 
-	log.Println("cuts")
-	for i := range cuts {
-		log.Printf("%.3f: %d\n", heights[i], len(cuts[i]))
+	log.Println("inside classes:")
+	for _, c := range classes {
+		log.Printf("\t%d\n", len(c))
+	}
+
+	log.Println("cuts:")
+	for i, c := range cuts {
+		log.Printf("\t%.3f: %d\n", heights[i], len(c))
 	}
 
 	// TODO: sew segments together.
@@ -130,9 +151,12 @@
 
 	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))
 
@@ -142,14 +166,14 @@
 		min.X, min.Y, min.Z,
 		max.X, max.Y, max.Z)
 
-	heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z)
-
-	log.Printf("classes: %d\n", len(heights))
-
+	start = time.Now()
 	tri, err := octree.Triangulate(xyz)
 	check(err)
+	log.Printf("triangulation took %v\n", time.Since(time.Now()))
 
-	start := time.Now()
+	log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)
+
+	start = time.Now()
 	intersectTriangles(tri, heights)
 	log.Printf("intersecting triangles took %v.\n", time.Since(start))
 }