# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1569496136 -7200
# Node ID e01b6e7cbe7d470c8e915175abbcf150f9feb57e
# Parent  99928231b679dc3f55f9433e3469280e8120b1b2
Detect isolated triangles as standalone areas.

diff -r 99928231b679 -r e01b6e7cbe7d cmd/isoareas/main.go
--- a/cmd/isoareas/main.go	Thu Sep 26 10:59:17 2019 +0200
+++ b/cmd/isoareas/main.go	Thu Sep 26 13:08:56 2019 +0200
@@ -18,6 +18,7 @@
 	"io"
 	"log"
 	"math"
+	"math/bits"
 	"os"
 	"strconv"
 	"strings"
@@ -83,7 +84,7 @@
 	}
 }
 
-func intersectTriangles(tri *octree.Triangulation, heights []float64) {
+func intersectTriangles(tri *octree.Triangulation, heights []float64) [][]octree.LineStringZ {
 
 	type indexedCut struct {
 		cut   octree.LineStringZ
@@ -97,18 +98,16 @@
 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]],
-			tri.Points[ti[1]],
-			tri.Points[ti[2]],
-		}
-		min := math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z))
+		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(t[0].Z, math.Max(t[1].Z, t[2].Z))
+		max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))
 
 		for j, h := range heights {
 
@@ -126,6 +125,7 @@
 			if min > h || max < h {
 				continue
 			}
+			t := octree.Triangle{v0, v1, v2}
 			cut := t.IntersectHorizontal(h)
 			if len(cut) >= 2 {
 				cuts[j] = append(cuts[j], indexedCut{cut, i / 3})
@@ -133,21 +133,39 @@
 		}
 	}
 
+	result := make([][]octree.LineStringZ, len(classes))
+
 	log.Println("inside classes:")
 	for i, c := range classes {
-		remain := make([]int32, 0, len(c))
+
+		var pb polygonBuilder
 
+		var isolated, inside int
+	allInClass:
 		for _, idx := range c {
-			// Eliminate triangles that do not contribute
-			// any border to the final outline.
-			if !inner(tri, idx, c) {
-				remain = append(remain, idx)
+			ns := neighbors(tri, idx)
+			w := where(ns, c)
+			switch bits.OnesCount8(w) {
+			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
 			}
 		}
-		log.Printf("\t%d (%d) %.2f%%\n",
-			len(remain), len(c),
-			100*float64(len(remain))/float64(len(c)))
-		classes[i] = remain
+
+		log.Printf("\t%d: inside: %d / isolated: %d\n",
+			i, inside, isolated)
+
+		result[i] = pb.polygons
 	}
 
 	log.Println("cuts:")
@@ -155,10 +173,21 @@
 		log.Printf("\t%.3f: %d\n", heights[i], len(c))
 	}
 
+	return result
+
 	// TODO: sew segments together.
 
 }
 
+type polygonBuilder struct {
+	polygons []octree.LineStringZ
+}
+
+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]
@@ -180,15 +209,29 @@
 	return false
 }
 
-func inner(t *octree.Triangulation, idx int32, indices []int32) bool {
-
-	for _, n := range neighbors(t, idx) {
+func where(neighbors, indices []int32) byte {
+	var mask byte
+	for i, n := range neighbors {
 		if n < 0 || !contains(n/3, indices) {
-			return false
+			mask |= 1 << i
 		}
 	}
+	return mask
+}
 
-	return true
+func dumpPolygonsSVG(
+	w io.Writer,
+	min, max octree.Vertex,
+	polygons [][]octree.LineStringZ,
+) error {
+
+	ratio := (max.X - min.X) / (max.Y - min.Y)
+
+	log.Printf("ratio: %.2f\n", ratio)
+
+	// TODO: Implement me!
+
+	return nil
 }
 
 func main() {
@@ -218,6 +261,8 @@
 	log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)
 
 	start = time.Now()
-	intersectTriangles(tri, heights)
+	polygons := intersectTriangles(tri, heights)
 	log.Printf("intersecting triangles took %v.\n", time.Since(start))
+
+	check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
 }